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ABSTRACT 


An imaging radar, like ISAR, offers a combatant the capability to perform 
long range surveillance with high quality imagery for positive target 
identification. Extending this attractive feature to the battle damage assessment 
problem (BDA) gives the operator instant viewing of the target’s behavior when 
it is hit. As a consequence, immediate and decisive action can be quickly taken 
(if required). However, the conventional Fourier processing adopted by most 
ISAR systems does not provide adequate time resolution to capture the target’s 
dynamic responses during the hit. As a result, the radar image becomes 
distorted. To improve the time resolution, time-frequency transform (TFT) 
methods of ISAR imaging have been proposed. Unlike traditional Fourier- 
based processing, TFT’s allows variable time resolution of the entire event that 
falls within the ISAR coherent integration period to be extracted as part of the 
imaging process. We have shown in this thesis that the use of linear Short 
Time-Frequency Transforms allows the translational response of the aircraft 
caused by a blast force to be clearly extracted. The TFT extracted images not 
only tell us how the aircraft responds to a blast effect but also provides 
additional information about the cause of image distortion in the traditional 
ISAR display. 
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I. INTRODUCTION 


A. ISAR FOR BATTLE DAMAGE ASSESSMENT (BDA) 

A radar targef image is generated from the reflectivity data collected by a 
radar platform as the target is observed from a set of viewing angles. For 
Synthetic Aperture Radar (or SAR), the radar platform moves to give the different 
target viewing angles, whereas Inverse Synthetic Aperture Radar (ISAR) makes 
use of the target’s angular rotations such as roll, pitch and yaw to collect the 
viewing angles for image generation. 

ISAR imaging opens up the new possibility of using radar measurements 
to construct high quality images for target recognition purposes while the target 
(usually aircraft and ships) is moving and the radar platform remains essentially 
stationary2. Operationally, ISAR offers enhanced target recognition at ranges 
that are often not achievable by electro-optic systems or Forward Looking 
Infrared (FLIR) alone. An extension of this long-range target recognition 
capability of the ISAR, apart from deploying for surveillance purpose, is to use it 
to perform Battle Damage Assessment (BDA). 

The significant advantage of using ISAR for BDA (over other methods) is 
its ability to provide the operator instantaneous viewing, in highly resolved detail, 
of the behavior of the target as a kinetic energy weapon hits it. This gives the 
operator immediate evaluation of the extent of the damage caused by the hit. 
Follow-up counter actions can be quickly taken if needed (all being performed 
beyond the target’s weapon range). 


B. RADAR TARGET IMAGING USING ISAR 

To construct a radar image from the reflectivity data, most ISAR systems 
process the data by performing a two dimensional (2D) Fast Fourier Transform 


"I The target can be any object of interest to the radar operator. 

2 The platform can be in motion but its motion has to be accurately measured (usually by the 
combination of Inertial Navigation System (INS) and Global Positioning System (GPS), and 
compensated accordingly. 
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(FFT); the first FFT generates the synthetic range profile of the target, and the 
second FFT extracts the Doppler information for each range cell. Ideally, the 
resultant 2D Range-Doppler plot following this process would closely resemble 
an outline of the target, arising from the various scattering points around the 
target. The quality of the generated ISAR image is directly related to the 
resolution imposed by the processing system. Fligh range resolution can be 
achieved through pulse compression, and Stepped Frequency Waveforms 
(SFW) are the most common methods used today. Fligh cross-range resolution 
depends strongly on the variation in viewing angles (equivalently Doppler 
resolution) which, in turn, increases with higher angular rate of target rotation and 
target dwell time. 

Flowever, a subtlety that affects the quality of the ISAR image results from 
the requirement of minimum variation in the range over the different viewing 
angles between the radar platform and the target for cross-range synthesis. 
Significant variation in range (due to the presence of translational velocity 
components, non-linear or large variations in angular velocity, or combination of 
both) will introduce an additional quadratic phase error that tends to distort or blur 
(known also as the defocusing effect) the image during the signal processing 
stage. This phenomenon is called “Range Walk” by the SAR imaging 
community, and it acts to impose an upper limit on the effective observed target 
rotation interval and—because this interval is determined by rotation rate—range 
walk also bounds the integration time. 

Normally, if the radial velocity remains essentially constant, the motion of 
the target can be accurately estimated (and compensated for) to reduce the 
effects of range walk. The angular velocity, on the other hand, is needed for 
ISAR imaging and is generally assumed to be constant and reasonably small 
over the entire data integration time. In most (aircraft or ship) target’s roll, pitch 
and yaw behaviors, this is approximately valid for the observation period. The 
converse is not true though if the target experiences any acceleration. 
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C. RESPONSE OF TARGET BEING HIT 

While the target is being hit, the blast force from the explosion of the 
weapon—or weapon impact itself—acting on the target structure will cause it to 
undergo translational and rotational motions, as well as structural vibration. The 
time duration of this force acting on the target is often very short (on order of 
milliseconds) but the impulse is large enough to cause significant alterations in 
both radial and angular velocity of the target. As a consequence of this short 
impulse period, the radar may not have adequate time to compensate for these 
changes in velocities and hence errors are introduced into the image processing. 
These errors could account for the distortion in most ISAR image observed when 
the target is being hit. 

In addition to the translational, rotation and vibration effects, the sharp and 
high velocity fragments from the weapon as it explodes, or possible a Kinetic 
Energy (KE) round itself, are likely to penetrate the target causing it to break up. 
The smaller broken target pieces now begin to scatter and reflect 
electromagnetic signals back to the radar as they move at much higher velocities 
(from the energy acquired) than the target (which has larger mass). The returned 
scattered EM signal from these pieces can behave like point scatterers, dihedral, 
trihedral or multi-facet corner reflectors, or specular reflectors depending on the 
dynamic nature of the break-up process. It is difficult to predict target structure 
changes without an in-depth analysis based on explosive and structural impact 
studies. In most cases, numerical simulations such as finite element analysis are 
often used in these studies and the results obtained are only specific to the 
particular scenario defined. Due to time constraints, this topic will not be pursued 
in the thesis work. 

D. JOINT TIME-FREQUECY METHOD OF ISAR IMAGING FOR BDA 

Time-frequency methods for signal analysis (also known as the Time- 
Frequency Transform or TFT) are similar to conventional Fourier Transform (FT) 
schemes except that the FT is taken after the auto-correlation of the signal with a 
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weighted window function. The advantage of this analysis method is that it allows 
decomposition of the frequency spectrum of time-varying signals in shorter 
window time frames, and provides an added dimension for examining the 
dynamic behavior of the signal as it varies over time. The size and shape of the 
weighting window function can be altered to fit the specific analysis requirement 
for the signal. 

Time-frequency methods in signal analysis have aroused much interest in 
the signal processing community. The initial application was in spectrogram 
analysis of complex human speech. However, in recent years, further extension 
of this time-frequency method to ISAR was actively advocated by Chen [2] and 
Jae and Thomas and Flores [3] — particularly in the area of image generation, 
motion compensation and micro-Doppler target vibration studies. In image 
generation, Chen [2] has reportedly used the TFT to construct ISAR image of 
targets with high rotation rate by helping to overcome the migration of individual 
scattering points from one range cell to another. Both Jae, Thomas and Flores 
[3] and Chen [2] have also proposed the use of TFTs to resolve multiple targets 
appearing in ISAR images. 

So far, the possibility of using the TFT to examine the dynamic behavior of 
a target in an ISAR image as the target is being acted upon by an impulse blast 
force has yet to be explored by other TFT-centered research efforts. If the target 
is assumed to be a rigid body, it would experience sudden translational and 
rotational motion variation under blast wave loading. These motions will be 
measured during the data collection process if the radar is continuously staring at 
the target. This thesis aims to apply the appropriate TFT window function to the 
collected ISAR data and demonstrate the possibility of extracting these target 
motions using TFT methods of image construction. 
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E. THESIS OUTLINE 

The thesis is organized in the following manner: 

Chapter II serves to develop the theory of ISAR image processing. This 
chapter will touch on how high quality image resolution is achieved, and the 
conditions necessary to maintain the quality of the image generated. 

In Chapter III, the effects of blast waves interacting with a target are 
described with the goal of modeling and simulation. The subject target for this 
study is an aircraft which is treated as a rigid body to simplify the analysis. The 
translational and rotational motions resulting from force interaction are examined. 

Chapter IV introduces the Time-Frequency Transform (TFT) Method of 
analysis, and how it can be applied to ISAR signal processing and image 
generation. The discussion in this chapter is largely based on the work of Chen 
[2] in TFT method for ISAR imaging and Jae and Thomas and Flores [3] for 
motion compensation. 

To demonstrate the feasibility of using TFT methods for ISAR imaging, 
and their advantage over conventional FFT based methods for BDA (oweing to 
the dynamic behavior of the aircraft as it is hit), a simple simulation program is 
developed. Chapter V describes how this program is put together for modeling 
the radar target returns from an aircraft, the effect of the blast loading on the 
aircraft and the construction of ISAR images using TFT concepts. This chapter 
also discusses the parameters and conditions used in the simulation. 

Chapter VI follows on to discuss the results of the simulation and its 
significance. 

Chapter VII is the concluding chapter and summarizes the major findings 
of the thesis work and provides recommendations for similar future work. 
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II. INVERSE SYNTHETIC APERTURE RADAR (ISAR) 
THEORY AND SIGNAL PROCESSING 

A. RADAR IMAGING AND ITS ASSOCIATED RESOLUTIONS 

Radar has long been recognized for its ability to “detect, locate and 
identify targets at great distance and in all kind of weather conditions” [1] and, to 
date, it is hardly matched by other sensors. Compared to Infrared (IR), visible or 
Ultra-violet sensors, the longer wavelengths (by four to six orders of magnitude) 
allow the radar energy to propagate through most of the atmosphere with little 
loss, even in poor weather conditions. However, the long wavelengths limit the 
radar’s imaging resolution, which is proportional to the ratio of the wavelength to 
the aperture size, and hence affects the image quality. To achieve the fine 
image resolution in optics, a lens system with aperture size of up to several 
centimeters is adequate. But a radar system would require an antenna aperture 
diameter to stretch tens of kilometers to attain the same resolution. Thus, 
alternate methods to accomplish this must be found. 

Inverse Synthetic Aperture Radar (ISAR) combines the use of pulse 
compression, high Pulse Repetition Frequency (PRF) and the target motions 
(particularly its angular rotation) to generate high range and cross-range 
resolution target images. Range resolution is associated with the width of the 
target pulse in the time domain. Hence, narrower pulses will have better 
resolution. But such radar systems require higher peak power to ensure enough 
energy returned from the target for signal processing. To mitigate this, pulse 
compression is mostly employed and for ISAR a common technique uses 
Stepped Frequency Waveforms (SFW). High cross-range resolution can be 
achieved by observing the phase variations of the signals returned from different 
target viewing angles over a coherent processing or integration period. It can be 
shown, under rather stringent conditions (in subsequent sections), that the 
changes in phase due the rotation of the target can be mapped in one-to-one 
correspondence to the cross-range position of the target scattering points. 
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B. ISAR SYSTEM ARCHITECTURE 

Figure 1 below shows the architecture of a simple SFW ISAR system. At 
the RF stage, a series of high PRF narrow pulses is used to modulate the RF 
carrier that varies in frequency steps of Af at every pulse interval, T 2 , which is 
reciprocal to the PRF. A burst of transmitted pulses consists of N stepped 
frequency pulses varying from fc to fc + (n-l)Af, where n = 1, 2, N.. For ISAR 

image processing, M bursts are required to form an image. The size of the 
image is thus determined by \he M x N sequence of the pulses, and the 
integration time, T, needed to construct the image (if we ignored the processing 
delays) is then T = M xN x T 2 . The M xN SFW pulse sequence represents the 
Doppler-Range mapping for the image. 



f] - Intermediate frequency 
fc - Carrier Frequency 
Af - Stepped frequency 
1,2, ..., N 


(a) RF transmission and detection stages (After: Figure 5.2(b), pp 202 , Wehner [4]) 
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Complex I and Q samples 



M Time Histories 


(b) Signal and image processing stages (After: Figure 5.2, pp 100 , Chen [2]) 


Figure 1 A simple illustration of SWF ISAR architecture. 

The returned echo signal captured by the receiver immediately after each 
pulse transmission is down-converted to an intermediate frequency (IF) and 
passed along for quadrature detection. The quadrature detector detects the 
baseband signal in the form of in-phase (I) and quadrature phase (Q) signals. 
These signals are then sampled and held as an M x N data array for signal 
processing. The sampled output data represents the time-history of the target’s 
reflectivity in the frequency domain. 

At the signal processing stage shown in Figure 1(b), corrections for radar 
system “phase and amplitude ripples” (Wehner [4]) are first applied, followed by 
range alignment, translational and rotational motion compensations in the form of 
range and Doppler tracking. Motion compensations are employed to minimize 
the induced quadratic phase error and hence “Range Walk” effect. The actual 
image construction begins with the synthetic range profile processing which is 
simply formed by taking an inverse Fourier transform of the N pulses for each m 
burst where m = 1, 2, M. Cross-range information is then extracted by another 
Fourier process - a direct Fourier transform this time - along the column of M 
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bursts for each n range cell index. In practice, the synthetic range profile and 
cross-range (Doppler history) information are simultaneously processed by taking 
a two-dimensional Fast Fourier Transform^ (FFT) of the compensated data 
array to lessen the computation load. 

C. THEORY OF MOVING TARGET IMAGING IN ISAR 

1. ISAR Geometry and its Relationship with the Returned Signai 

For illustration purposes, consider a simplified two-dimensional (2D) 
geometry between the radar and an airborne target as shown below in Figure 2. 
The aircraft is assumed to have translational and rotational motion only in the 
two-dimensional plane relative to a stationary radar platform. It can be easily 
extended to the actual three-dimensional case to include the effect of the motions 
along and about the other axes with some careful considerations. 


(o rotation 



y cos 6* 


Figure 2 ISAR geometry in two-dimensional space. 

(After: Figure 5.1, pp95, Chen [2]) 

3 It can be shown that the process of IDFT and subsequently DFT is equivalent to taking a 
two-dimensional FFT with the exception of a phase factor in front. In most radar signal 
processing, we can ignore this phase factor since our interest is largely on the power spectral 
density or absolute intensity of the returns. 
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Following Chen [2], we designate the radar to be located at the origin of 
the U-V coordinate system, known as radar coordinates. The target can be 
described by a collection of scattering points (or, “scattering centers”) in its own 
reference frame, x-y, known as the target coordinates. The center of rotation of 
the target is located at the origin of this x-j coordinate system. To describe the 
rotation of the target independent of the target scattering center distribution in x-_y 
coordinates, another reference coordinate system denoted as X-Y is introduced 
with its origin coincident with the target coordinates. The two axes in X-Y are 
parallel to U and F axes respectively so that a pure translational transformation 
can be used to describe their relationship. Then, it follows that the 
transformation of x-j reference to X-Y coordinate can be described by a pure 
rotation. Also, to simplify the illustration further, we assume a = 0 (the radar 
target line of sight (LOS) is along the U-axis of the radar reference). 

Assume that at t = 0 the target range (the distance from the radar to the 
center of rotation of the target) is R. Hence, the distance between the radar and a 
scattering point, P(x,y), can be expressed as 




(i? + xcos^o -ysmOf^'f +(jcos6*o + xsin6*o)^ 


1/2 


= +x^ +_y^ + 2i?(xcos6’o -_ysin6*o)j 

= i? + xcos6’o - jsin6*Q 


( 2 . 1 ) 


where 6o is the initial target rotation angle. 

If the target possesses both translational and rotational motion described 
by an initial angular rotation rate [co), initial radial velocity [Ft) and its higher order 
initial angular {}) and radial acceleration [at) and so on, the target range and the 
rotation angle as a function of time are then 

R{t) = R^ + F,t + ^a,t^ +.... ( 2 . 2 ) 
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(2.3) 


1 9 

0(t) = + (Ot + —]^t +_ 

Thus, at time t, (2.1) becomes, 

Rp = R(t)+ x cos 0(t)-y sin 0(t) (2.4) 

For simplicity, if we assume the reflectivity function of target, p(x,y), to be a 
pure function of the target geometry and that it does not change with the above 
parameters, then the base-band of the returned signal from the scattering point, 
P(x,y), will consist simply of the product of the reflectivity function and the phase 
variation due to the distance, Rp\ 

Spit) = pix,y)cx^\-j27tf^ ^ p{x,y)cx^{-jy/iRp,)] (2.5) 


2. Signal Representation 

We next look at the signal representation of the returned echoes for the 
imaging process. At any given time t, the radar receives a collection of returned 
signals from all the scattering centers of the target. The received base-band 
signal can be described as a superposition of these signals: 

-j In X dxdy (2.6) 

Substituting (2.4) into (2.6), we obtain 

Spit)^e^ \ \ pix,y)cxy)\-j2n[xf^it) + yfyit)~^dxdy (2.7) 

-00 -00 

where 

f^{t)=^coseit) and 4(0 = -^sin^(0 
c c 

Apart from the pre-factor (due to the translational motion) in (2.7), the 
base-band received signal at a given time can be seen to be a two-dimensional 
Fourier transform of the target reflectivity function. This is the basis for the 
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((0= j j /?(x,j)exp< 



method of extracting the target reflectivity information by applying a direct two- 
dimensional inverse Fourier transform or FFT on the collected data. 


3. ISAR Data Collection and Processing 

In almost all ISAR systems, the signal processing is performed discretely. 
To gain more intuition into this procedure, it is convenient to consider the discrete 
representation of the base-band signal rather than its continuous form. For our 
earlier discussion on ISAR architecture, we set the frequency dependence for 
each m burst of pulses to increase from fc in steps of 4/’from pulse to pulse. 

fn ~ fc +(^“ 1 ) 4 /^ where n = 1, 2, N (2.8) 


Also, the base-band signal will be time-sampled at the interval corresponding to 
the pulse interval (T^) plus the range of the target center. Accordingly (Jae, 
Thomas and Flores [3]), the discrete sampled frequency and time can be 
represented as 

2i? 

t^i^=\n + {m-\)N^T 2 ^ -^ where m =1, 2, M (2.9) 


Therefore, (2.2), (2.3) and (2.7) can be expressed in terms of (2.8) and 
(2.9) as 


Sj^(m,n) =e " p(x, y) exp {- 72 ;z- [x/^ (m, n) + yf^ (m, n)]} dxdy 


lR{m,n) 00 oo 


R(m,n) = R,. + V,t„ „ + f -i-.... 

\ ^ r () t m,n 2 t m,n 


2 /„ 

c 

2 /„ 

c 


( 2 . 10 ) 




fM^n) = -^sin6'(^o +••••) 
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Figure 3 Row-column decomposition for synthetic ISAR imaging. 


Figure 3 shows the image-processing sequence for the M x ISAR 
complex-valued sampled data matrix array. The image construction process 
begins with synthetic range profile extraction by an Inverse Discrete Fourier 
Transform (IDFT). This operation on each m and n point in the Mx matrix is 
given as such (Jae, Thomas and Flores [3]). 


h(m,n) = —X s^(m,i)e 
^ 1 = 0 


.2 Trni 

t - 

N 


where n = 1, 2, N 


( 2 . 11 ) 


The resultant M x N matrix is analyzed by another Fourier process - the 
Discrete Fourier Transform (DFT) - to gather the cross-range information. If we 
denote the output of the transformed point in the resultant MxA^ matrix by D(m,n), 
we obtain 


1 M -1 

D(m,n) = —Y h(k, 
M 


n)e 


. lirkm 


where m = 1, 2, M 


( 2 . 12 ) 
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The magnitude of the elements of resultant complex matrix forms the ISAR 
image. 

D. ACHIEVING HIGH RANGE RESOLUTION THROUGH STEPPED 

FREQUENCY WAVEFORM (SFW) SIGNAL PROCESSING 

Most radars today attain high range resolution through various forms of 
pulse compression. For a pulse compression radar, the achievable range 
resolution is inversely proportional to the system bandwidth, and thus, to gain 
high resolution the entire radar system — from transmitter through receiver and 
including data collection — must possess the large bandwidth associated with 
the desired resolution. Designing a radar to have such wide bandwidth at all its 
stages can be an expensive and a rather complex issue. 

SFW avoids this stringent requirement by simply varying the carrier 
frequency prior to transmission in N frequency steps associated with the desired 
range resolution. The same stepped carrier frequency is also used as the 
reference to the mixer to down-convert the returned signal back to the IF stage. 
Consequently, the IF stage, the detection stage and later processing stages no 
longer require large bandwidth to maintain the resolution. In addition to relaxing 
the wide bandwidth requirement for the other radar system stages, the SFW 
technique also has the advantage that the returned base-band signal is a 
spectral representation of the target’s reflectivity with range-invariant magnitude. 
This is a desirable property for range tracking (Wehner [4]). 

In doing so, on the other hand, the returned signal no longer preserves its 
direct phase relationship with the range profile of the target. To obtain the target 
range profile an inverse Fourier transform is needed to get back the target range 
profile. This process is known as synthetic range-profile processing. 

1. Synthetic Range Profile Processing 

Figure 4 shows an SFW transmission and its associated returned echo 
pulse. From Wehner [4], each received echo pulse is equivalent to a replica of 
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the transmitted SFW with a time delay of T(t) = 2 (Rp - vtt) / c. After passing 
through RF, IF, and the quadrature detection stage, the carrier component is 
removed leaving only the base-band signal with a phase output of Wn(t) = 27fnT(t). 
This signal is, in turn, sampled at = nT 2 + Tr + 2Rp/c where is the receiving 
system transfer delay, Rp represents the range of a point target (it is often 
estimated through range alignment), and c is the velocity of wave propagation 
(Wehner [4]). 


Transmitted Pulses 



fc+(i-mf 



Figure 4 SFW Pulse and echo pulses. 
(After: Figure 5.5, pp 204, Wehner [4]) 


The sampled output from both the I- and Q- channels can now be written 
in complex form where An is the signal amplitude (which is usually normalized for 
ease of illustration (Wehner [4])): 

(2.13) 


and 
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Wn = -2^/„ 


(2.14) 


2v 

f 2i?^ Y 

p t 

nT^+T,.+ 

c c 

y ^ ). 


The N sampled output data for each burst are then Fourier transformed by 
IDFT or FFT to obtain a series of synthetic complex range-profiles of the target, 
Hi, where / is the slant-range position. Following Wehner [4], the IDFT is 
expressed as 

N-\ 

, 0 < / < iV -1 (2.15) 

1=0 


Putting (2.13) and (2.14) into (2.15), and assuming zero target velocity, the 
normalized synthetic response of the target is 


u V • r 2^ ,. ^ ^ 2i? 


i=0 


N 


(2.16) 


Replacing^; =fc + iAf in (2.16), this becomes 


Hj = exp 




2R^'\N-X 




1=0 


Ini 


N 


-2NRAf 


+ / 


(2.17) 


This equation can be further simplified using the same analysis as for antenna 
array theory (from Wehner [4]) to become 


Hi = exp 




2R. 


sin Try 


2 sin 




exp 


J 


.N - \ Ijiy 


2 H 


(2.18) 


H 


where 


-2iW?^A/ ^ , 

_y =-+ / 

c 


(2.19) 


It is noted for (2.18) that the resultant response represents the range point 
spread function (PSF) of the target in the range cell index / where the highest 
peak indicates the range position of the target point. The collection of all the 
range cell indices from / = 0, I,..., H - l\s known to be the synthetic range profile 
of the target. 
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2. Unambiguous Range and Range Resolution 

From both (2.18) and (2.19), the peak response occurs when y = 0, +/-N, 
+/-2N, +/-3N,... corresponding to the range given by (Wehner [4]) 

2NAf 2NAf 2NAf 


The unambiguous range length is then: 

AD c(L + N) cl cil + 2N) c{l + N) c 

AR = — - - — or — - - - or... =- 

2NAf 2NAf 2NAf 2NAf 2Af 


( 2 . 21 ) 


and the range resolution (using Rayleigh’s resolution distance criteria^) 
corresponds to the width of the PSF of each target point, which is 


AR 


p 


c 

2NAf 


( 2 . 22 ) 


Both parameters depend on the frequency step size but the range resolution can 
be enhanced by placing a higher number of pulses within each burst. 


3. Effect of Target Velocity 

In the above, the synthetic range response of the target is obtained for a 
target with zero velocity. Flowever, if target linear velocity is present, the 
synthetic range response becomes 




2R. 



f 

^^2 + + 


V 



(2.23) 


In this case, the range resolution can no longer be simply expressed as in (2.22). 
Moreover, analytically deriving the range resolution ARp from (2.23) is also not a 
trivial problem. Through numerical solution, Wehner [4] has reportedly shown 


^ Rayleigh criterion defines the resolution between two point spread functions to be the 
overlap of one highest maximum peak with the first minimum of the other highest peak. For a 
sine function, this occurs at the first null in which the argument in the sine function is n. 
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that the higher velocity can cause the width of the PSF to broaden, thus resulting 
in poorer range resolution. 


4. Range Offset and Range Walk 

The presence of nonzero translational target velocity has the 
consequence of creating range offset and range walk effects, and therefore 
distorts the radar image. When ISAR constructs a radar image from the collected 
data, it presumes that the range of every target scattering point does not vary 
during the integration time interval. Range offset and range walk result when this 
presumed condition fails to be true due to the presence of radial velocity 
components. The extent of the range walk and range offset can be quantified by 
the number of cells {AN) that the range is shifted over the entire time integration 
period (7) for the radar time delay resolution oi At = 1 /NAf. Thus, from Wehner 
[4], AN \s expressed as 

AN = — = NAf St (2.24) 

At 


where St represents the time delay shift of the target due to the presence of 
radial velocity over the period where St is evaluated. This delay shift is 
determined from the total group delay and can be derived by taking the derivative 
of (2.14) and subtracting the time delay 2Rp/c. In general, ^ris cumulative over 
the number of bursts, m, and is given by Wehner [4] as 


St = 


2v, 


cAf 


f 


fcT2+¥ 


2i? 

r -\ - \-Nmf 




(2.25) 


Putting (2.25) into (2.24), and noting that NmT 2 » Tr + 2R / c is true in most 
cases, we obtain 

AN \ f+ NmAf) V, (2.26) 

For the entire image frame time, we let m=M, and obtain 
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AN- 


(2.27) 


2(NT,) ^^ ^ 2M(NT,)NAf 


The first term in (2.27) represents the range offset which increases with 
higher values of N, higher carrier frequency and lower PRF. The second term is 
the result known as range walk and is cumulative from profile to profile (burst to 
burst) depending on size of M. Unlike range offset, the non-constant range shift 
due to the cumulative nature of range walk tends to move the same scattering 
points out of the range-cell column (known as range cell migration) and affects 
the cross-range information extraction process (which, in effect, takes DFTs of 
each range-cell column). This will generally result in distortion of the ISAR image. 
The choice of Af also contributes to the range walk effect. For better range 
resolution, we require high 4/" but range walk will tend to be amplified by higher 
Af- 

To correct the distortion caused by range cell migration of the scattering 
points, the phase error due to the target velocity has to be mitigated. The 
simplest method for eliminating this error is to multiply the sampled output data 
by an equivalent complex factor (Wehner [4]) 


g = exp 




2v. 


^ 9 5 

V ^ 


(2.28) 


where v, and are estimated from the Doppler and range tracking of the 

sampled output data. This is known as translational motion compensation. 
However, it is evident from (2.28) that the effectiveness of this kind of 
compensation relies heavily on the precision of Doppler and range tracking 
techniques in order to get an accurate estimated of the target range and its radial 
velocity. Compensation methods get more complex if the target velocity varies 
significantly within the integration period. 

There are many methods today that can yield accurate motion estimation. 
One of the known methods is the entropy measurement but this method is 
computationally intensive (Jae and Thomas and Flores [3]). Other techniques 
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include the Cramer-Rao Burst Derivatives Approach (Jae and Thomas and 
Flores [3]), the Phase Difference method (Wehner [4]), the least square motion 
parameter estimation (Jae and Thomas and Flores [3]), and even time-frequency 
analysis approaches (Jae and Thomas and Flores [3] and Chen [2]). 


E. ACHIEVING HIGH CROSS-RANGE RESOLUTION USING SYNTHETIC 

APERTURE 

1. ISAR Theory from an Aperture Viewpoint 

The method by which ISAR achieves high cross-range resolution is not 
much different than that of conventional SAR and, in both cases, the bounds are 
linked closely to the angular extent over which the radar views the target. For an 
unfocussed SAR, it is the maximum linear path length before the quadratic phase 
error becomes significant to create the defocusing effect and degradation of 
resolution. (This phase error results from the difference between the linear path 
and the curvature of the range to a point target.) In the case of ISAR, the same 
behavior holds except that the path length is now a function of the target angular 
rotation rate and the dwell time (which is also usually the image coherent 
processing or integration time). 

Consider a stationary radar viewing a rotating target (see Figure 5(a)) over 
an angle segment W (aperture) over the integration period (we assume that the 
translational motion is perfectly compensated). This geometry can be viewed as 
equivalent to the radar (like SAR) in motion observing a non-rotating target as it 
passed the same angle segment as shown in Figure 5(b). We assume that the 
radar is tracking the target in the far-field so that the angular extent of the 
aperture, projected as path length, L, is L « R. Also, the length of the observed 
target, y, is generally _y << L. In such cases, from Figure 5(c), the two-way phase 
advance of the target scattering point, cp, at azimuth position, from the radar 
boresight is 

Ajr 

(p{x)'^lkdr = —x%m.(j)^ (2.29) 
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Stationary Radar 



(a) ISAR (After: Figure 7.3(a), pp 344, Wehner [4]) 



(b) SAR equivalent (After: Figure 7.3(b), pp 344, Wehner [4]) 



(C) Approximation for L«R and y«L (After: Figure 7.3(c), pp 344, Wehner [4]) 

Figure 5 SAR and ISAR equivalence. 
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The term sin (j)y approaches y / R\ox y « R, and the projected length x 
approaches R^ under small angle approximation. Hence, (2.29) can be 
expressed in terms of ^/iand becomes (p(4>) = 4n <f)y / X. Using similar arguments 
as in linear antenna array theory, the normalized integrated response over the 
integration angle, % behaves like a point spread function (PSF) about the cross¬ 
range position, y: 

= = (2.30) 

J-T ^y/y 

Using the Rayleigh resolution criteria, we can determine the cross-range 
resolution between two point targets to be Arc = X / 2^. For small integration 
angle ISAR, we set *F= coT, where denotes the target angular rotation rate and 
ris the integration time. Then, the cross-range resolution is given by 

Ar^=— (2.31) 

" 2a)T 


Like SAR, the size of the target imposes a maximum allowable 'F for 
which this cross-range resolution can be achieved. It can be shown (see Wehner 
[4]) by setting the maximum allowable phase derivation to be no more than n/8 
rad as the criteria for maintaining focus, that the maximum integration angle 
before defocusing occurs is 


¥ 


max 



(2.32) 


where the maximum target length is L, and to maintain focus it is necessary that 


L< 


2X 


(2.33) 


2. Range Doppler Imaging and Cross Resolution 

The cross-range resolution in the above situation is obtained by drawing 
an analogy to SAR in terms of geometry. The aim is to emphasize the 
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importance of maintaining a small integration angle, *F, within the linear regime 
for which the quadratic phase is not significant. Unlike SAR, however, in ISAR 
this angle is affected by the target rotation rate and is often not a known priori — 
nor can it be generally measured directly by the radar. 

In most — if not all — ISAR systems, Doppler processing techniques are 
directly applied to the signal to obtain the cross-range information. This 
technique avoids the need for exact knowledge of the target rotation rate, co, as 
long as the product of coT remains within the small integration angle limit. The 
relationship between target rotation, the point scatterer’s cross-range position 
and its resulting Doppler shift can be understood by referring to Figure 6. If we 
ignore the Doppler component due to translational motion (it will be shown later 
that this component is just a linearly added term to the Doppler shift), then the 
Doppler shift due to the instantaneous velocity of the scattering point at a 
distance from the center of target, arc, is 



Figure 6 Radial velocity from a point scatterer on a rotating target. 
(After: Figure 7.6, pp 351, Wehner [4]) 


From Equation (2.34) (after some manipulation), the cross-range resolution is 
seen to be directly related to the Doppler resolution of the radar, Afu. 
Furthermore, since Afu = 1/T, the cross-range resolution can be expressed as 


Ar 


-4/n =- 

r\ U r\ rri 

lo) 2a>T 


(2.35) 
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This expression is the same as that obtained earlier for the cross-range 
resolution. 

3. Doppler Frequency Shift Due to Target Motion 

The simple relationships above may not reveal much about the effects of 
target motion of the ISAR image. To give further insight, let’s consider the phase 
relationship derived earlier in Equation (2.5). An approximation of the Doppler 
shift induced by the target’s motion can be determined by taking the time- 
derivative of the phase term y/(Rpt) where Rp, ^ R{t) + xco?, 0{t) - y sm 0{t), 

R{t) = R^^ + V,t + \af+.... and 0{t) = 0Q+cot + \yt^ + . Neglecting the second 

order and higher terms of R(t) and 6(t), we obtain 
d 2 Z' 2 Z' 

=— i//(Rp^) = -l-^^{-xo[sin6*oCOS(yr + cos^oSin®r] 

dt c c 

-ya)\cos 6 q cos cot + sin sin or]} (2.36) 

fDTrans fDRot 

The overall Doppler shift can be seen as a summation of the Doppler shift 
due to translational motion,= 2fcVt/c, and that due to rotational movement, 
foRot, determined by the remaining terms. We can further expand sin cot = cot - coh^ 
/ 6 + ... and cos cot = 1 - coh^ / 2 + ... and ignore the higher order terms since the 
product of both co and t is generally very small, oV << 1 and ooi « 1. Thus, 
fuRot becomes 

2 f 

foRot =-^{-®[xsin6'o+jcos6’o]-®^r[xcos6’o-jsin6’o]} (2.37) 

The first term in (2.37) is used for ISAR imaging. But there exist quadratic 
phase terms even though the rotation rate remains constant. These quadratic 
terms vary with time and serve to introduce an additional phase error that can 
defocus and distort the image during its construction process. 
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4. Distortion Produced by Target Rotation 

Unlike errors due to translational motion, the phase error due to target 
rotation can cause cell migration in both range and cross-range dimensions and, 
hence, distortion in the ISAR image. From the 2D radar-target geometry given 
earlier, cell migration in the range dimension due to target rotation is most 
significant when the initial target angle is 6o = n/ 2. The number of cells moved in 
the range dimension over the integration period T is then (following Wehner [4]) 

AM' = aiA {i;,(0) - ^,(7-)} 

Rp (t) = + V,t + x 008(6*0 + J sin(6*o + cot) 


To account for time delay shift due to rotation alone, we set Vt = 0. 
Expanding the sine and cosine terms in Rp and setting 9o = n / 2, we can 
approximate AM’ as 


AM ' = — {y - xsin coT - y cos coT } 
Ar 


J_ 

Ar 


j-x 






^ {(oTf ^ 


xcoT 

Ar 


(2.39) 


Cross-range cell migration is most significant when 6o = 0. The number of 
cross-range cells by which the scattering point is offset from its correct position 
due to target rotation rate over the integration period T \s then 

AM=absMf„JO)-f,^,(T)} (2.40) 

A/ D 

Substituting (2.37) into (2.40) and letting 6o = 0 and Afp = 1 / T, \Ne get 

AM g ^ (2.41) 

c 

Earlier, we defined the cross-range resolution as Arc = A./2coT= c /2fcCoT, and so 
by substituting this result into the above equation, we obtain the same expression 
as in Equation (2.39). 
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(2.42) 




xcoT 

Ar 

C 


The extent of the cell migration is dependent on the target size, its rotation 
rate and the integration time. Also, higher resolution requirements over the same 
period T would mean increasing the number of range {N) and cross-range cell 
{M), and a corresponding increase in AM and AM’. This is similar to the effect of 
increasing the range resolution discussed in the previous section. 

Often, ISAR can be made to remain focused over the entire rotation angle 
T if the target phase drift in both range and cross-range (also known as blur 
radius) caused by target rotation is constrained to lie within one cell. From 
Equation (2.42), if we set AM < 1 for focused ISAR and for a>T = A / 2Arc, this 
yields 


x< 


2(Arf 

A 


(2.43) 


Equation (2.43) is almost the same as Equation (2.33) which was obtained 
by setting the maximal phase deviation to tt/ 8. The equation is determined by 
“the extremes of the target extent” (Wehner [4]), the required resolution and the 
frequency used. 
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III. BLAST LOADING ON AIRCRAFT 


A. BLAST IMPACT ON AIRBORNE TARGETS 

Explosive munitions, guided or unguided, are extensively used in most 
weapons today against all kinds of targets: aircrafts, ships, tanks, troops, 
infrastructures, etc. These munitions are attractive because of the enormous 
chemical energy that is released in a very short period when the explosive is 
detonated. The huge energy deposited into the medium creates a high velocity 
high pressure shock front that has damaging effects on the target structure — 
especially for buildings and fixed infrastructures. When used against airborne 
targets, the primary goal is to make use of this high velocity pressure wave to 
launch small pre-fabricated fragments and shrapnel which will penetrate different 
parts of the target body, causing a variety of damage mechanisms such as 
structural tearing, explosion of the fuel tank, loss of flight control and electronics, 
and injuring or killing the pilot. 

Owning to its finite inertial mass, the speed at which these fragments and 
shrapnel travel will in general be less than the shock-wave speed. As a 
consequence, the blast wave will interact with the target body before the 
fragments and shrapnel arrive. Most modern airborne targets are constructed 
from material with tensile strength able to sustain high aerodynamic and 
environmental drags of flight. Also, because of the degrees of freedom of the 
aircraft in air, the net effect of a blast force is more likely to cause translational 
and rotational motions (as if it is treated as a rigid body) than yielding its structure 
strength. A brief impulse moment can thus be captured by ISAR for Battle 
Damage Assessment (BDA) purposes. 

For the fragments and shrapnel that penetrate the target body, different 
damage mechanisms can occur depending on the nature of the penetration. If 
the fragments hit a fuel tank, a secondary explosion will occur and may cause 
further structural breakup. Alternatively, the fragments may cause the target to 
lose flight control by damaging the electronics, aero-structures, control 
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hydraulics, or even killing the pilot. In such cases, the target will display erratic 
maneuvers. Other scenarios include tearing out the tail and/or the wings, 
damaging the engine, detonating the weapons stores, etc. Although these 
effects are also captured by ISAR, they are in general difficult to completely 
account for even in computer modeling and simulations. 

In this thesis, we would like to briefly examine the effects of an explosive 
blast on the motions of the airborne target and demonstrate how it can be 
detected using ISAR. These imposed motions are highly dynamic and depend 
on strength of the explosion, the relative distance and direction in which the 
detonation occurs, and the surface area of the aircraft subjected to the blast 
loading. Very often, the solution to these complex motions can not be easily 
derived analytically. However, numerical solution of these complex equations of 
motion is possible. Hence, this chapter aims to provide discussions linking the 
theory of explosive blast in air, the loading of the blast force on the aircraft, and 
the aircraft’s translational and rotational responses through a set of complex 
equations of motion. In addition, a methodological approach is also proposed to 
create a computer model to solve these complex equations numerically. 


B. CHARACTERISTIC OF EXPLOSIVE BLAST IN AIR 

1. Formation of the Blast Wave 

Blast waves are formed as a result of the ambient atmosphere being 
forcibly pushed back by the high pressure gases produced from an explosion. In 
the initial formation, the shape of the pressure pulse from a conventional 
chemical explosion is highly dependent on the geometric construction of the 
explosive charge. However, as the pulses travel in the medium, the higher 
pressure portions will possess higher speed than the lower pressure parts 
(Kinney and Graham [10]). Progressively, the wave front becomes steeper as 
the pressure pulse move out and develops a discontinuity known as an explosive 
shock as illustrated in Figure 7. Under this condition, increasing amounts of 
energy will be needed to further alter the shock front shape. The front will largely 
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remain in this form and acts as an advancing blast wave system propagating 
away from the charge in the radial direction. 





(a) (b) (c) 

Figure 7 Development of explosive shock, (a) Initial pressure pulse, (b) successive 
configurations, and (c) Final form when fully developed. 

(After: Fig 6-1, pp 89, Kinney and Graham [10]) 


It is interesting to note that “the fully developed explosive blast-wave 
system is always formed with about the same general configuration no matter 
what is assumed for the initial positive pulse. That is, any initial configuration is 
lost and all blast waves at reasonable distances from the center of an explosion 
show similar wave forms” (Kinney and Graham [10]). 

The amount of energy released by the chemical reaction process from the 
explosion is an important factor in determine the intensity of the initial shock front 
and its subsequent impulse shape. In its initial form, the explosion process 
entails a vanishingly small volume of gases under high pressure. As the gases 
expand outward, the same amount of energy will need to overcome a larger 
volume of the atmosphere (which acts to resist the expansion process) hence 
lowering its resultant shock front intensity and ‘stretching’ the wave shape as 
shown in Figure 8. 
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Figure 8 Typical pressure-distance curves for successive times after an explosion. 

(After: Fig 6-2, pp 90, Kinney and Graham [10]) 


At some appreciable distance, the shape of the blast wave will depict a 
series of positive and negative phases of the pressure. In general, the first 
positive pressure phase is far more intense in causing damage than the negative 
ones and its subsequent rarefactions which are limited in magnitude to those of 
the atmosphere (Kinney and Graham [10]). 


2. Characteristic of the Biast Wave 



Figure 9 Typical pressure-time curves for explosive blast wave. 
(After: Fig 6-5, pp 99, Kinney and Graham [10]) 
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Figure 9 shows a typical pressure-time history characteristic curve of an 
explosive blast at a certain location away from the center of the explosion (after it 
is fully developed). At arrival time (4) after the explosion, there exists a sudden 
jump in the pressure intensity (known as peak overpressure) from its 
atmospheric level. An object at this location will be subjected to an 
instantaneous lateral force equal to the product of the overpressure and the 
projected area in the plane of the blast wave. However, this is not a stable 
condition, and the overpressure at the location will begin to decay exponentially 
as the volume of gases continues to expand outward, leading to smaller negative 
and subsequently smaller positive rarefaction phases. For most purposes, the 
main interest lies in the first positive phase of the blast wave where the damaging 
effect is considerably more significant than its later ones. 

There are four independent parameters that can be used to completely 
describe the first positive phase of blast wave shape: (a) the initial shock intensity 
specified by peak overpressure (po), (b) the duration of the positive phase of the 
blast {td), (c) the impulse (force-time product) per unit area [I/A), and (d) the 
arrival time (4) which merely introduces an overall time shift of the pressure-time 
history curve. 

An extension of the logarithmic decay expression from Kinney and 
Graham [10] linking these parameters to the blast wave shape can be given as 


p[t) = 


Po 

0 


^ t-t ^ 
V J 


exp 


a{t-tj 


where t^<t<tj 
otherwise 


(3.1) 


where p is the instantaneous overpressure at time t, and a is the waveform 
parameter related to the impulse per unit area. 
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3. Peak Overpressure 

Theoretical analysis of the peak shock overpressure uses the same 
approach as for a normal shock. It is, however, complicated by the spherical 
divergence of the shock front and its transient nature. Both theoretical studies (in 
some cases, through the use of computer solutions) and experimental 
measurements are usually correlated to obtain a reasonable expression relating 
peak intensity to the distance. Following Kinney and Graham [10], the 
empirically derived peak overpressure-to-distance relationship can be expressed 
as the ratio of the overpressure to the standard atmospheric pressure of 
1013.25mbars at 15deg C. 


Pa 


798 


1 + 




1 + 


0.048 




1 + 


I r z ^ 


0.32 


1 + 


1.35 


(3.2) 


where po/Pa is overpressure ratio and Z is the scaled distance of a charge weight 
at 1kg TNT equivalent. The scaling of Z to the actual distance will be discussed 
below. 


4. Arrival Time 

From Kinney and Graham [10], we noted that the velocity of the shock 
front is uniquely related to po/Pa which in turn depends on the distance between 
the observed point and the center of the explosion. Flence, given the peak 
overpressure and the distance over which it occurs, we can obtain the arrival 
time by noting 

= (3.3) 

dt 

where Ux is the shock front velocity, Mx is its Mach number, r is the radial 
distance, and Qx is the speed of sound in the undisturbed atmosphere. Solving 
(3.3) by separation of variables and integrating the Mach number from the center 


34 



of the charge (assumed to be rc = 0) to the distance the shock front arrived, we 
obtain 


1 f 1 , 

K = — 

0 


where the expression for shock velocity in term of Mach number (given by Kinney 
and Graham [10]) is 


M 


X 


ik + \)p^ 


where k = l A 


(3.5) 


5. Blast Duration 

The extent of the damage caused by the blast wave is also closely linked 
to the duration over which the force is applied. Since the impulse (equivalent to 
the area under the blast wave shape) of the positive phase is much higher than 
its subsequent phases, it is often used as an index for the blast duration although 
in “some cases the negative phase duration can be twice as long” (Kinney and 
Graham [10]). Thus, the blast duration at an observed point is defined to be “the 
time between the passing of the shock front and the end of the positive pressure 
phase as marked by a zero overpressure” (Kinney and Graham [10]). 

The zero overpressure condition is a characteristic of the sound wave 
since the shock front velocity will decrease to the speed of sound in the medium. 
Hence the blast duration is also dependent on the shock velocity. Using a similar 
approach to that for deriving the arrival time, Kinney and Graham [10] give the 
relationship between the blast duration and the distance to which the shock front 
begins to develop as 
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where the term td/ is the duration in msec for a 1 kg TNT explosion and Z is 
the same scaled distance as in the peak overpressure expression. The actual 
blast duration can be obtained by multiplying (3.6) by where W is the weight 
of the explosive used in kg. 

6. Impulse Per Unit Area 

Impulse has “the dimensions of the force-time product” and can be 
obtained graphically from the area under the first positive phase of the blast wave 
characteristic. This factor determines the extent of damage to the target 
structure. From Figure 9, we can see that — apart from the peak overpressure 
and blast duration — the rate of decay of the overpressure (known as the 
waveform parameter, a) also affects the amount of impulse force acting on the 
target. A rapid decay blast wave with the same intensity will have less damaging 
effect on the target than the slow decaying ones. 

Flowever, it is often difficult to determine a from analytical analysis or 
direct measurement. In contrast, the impulse per unit area {I/A) can be easily 
obtained by experiment, and the empirical relationship between I/A and the 
scaled distance (from Kinney and Graham [10]) can be expressed as 



With Equation (3.7), the waveform parameter in turn can be determined by 
rearranging the following integral expression 

IIA=']‘ pdt = pj\--^(l-e‘-'^ (3.8) 
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7. Hopkinson Scaling Law 

Scaling of parameters allows the quantitative characteristic of any 
explosion (especially those drawn from experimental data) to be used for 
solutions to more general explosive blast wave problems. The scaling law is 
fundamentally based on the geometrical similarity between the reference and 
actual object. The general principle for explosive scaling is to consider the two 
objects to be spherical. Then, by geometry, the ratio of their volumes is 
proportional to the third power of the ratio of their diameters, and if the two 
objects have uniform density distribution, their mass ratio also follows the same 
third power rule. 

By extending the above argument, the blast wave of two explosions can 
be said to be identical if the ratio of their distance is proportional to the cube root 
of their energy release. If we take the medium into account, it has been shown 
by Kinney and Graham [10] that the scaling for distance can be expressed as 
follows: 


{actual _ dis tan ce){atmospheric _ density)^'^ 
{energy _ release)^'^ 


{scaled _ dis tan ce) 


(3.9) 


Since the energy released is related to the velocity of the shock front, 
which in turn depends on the peak overpressure ratio and hence the atmospheric 
density, the above expression can be simplified to become 

^ L X {actual _dis tan ce) 

Z = {scaled _ dis tan ce) = — -- (3.10) 


where fd is known as the atmospheric transmission factor and accounts for the 
difference between the actual and standard reference atmosphere (designated 
by subscript o). /^can be found using 


fd = 


\Po J 


1/3 




1/3 


yPoj 


f T 


J 


(3.11) 
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The overpressure is indirectly determined by the scaled distance relation, 
and the actual overpressure is simply 

Overpressure = overpressure _ ratio x atmospheric _ pressure (3.12) 

This is the situation of the Hopkinson Scaling Law (Ref from Kinney and Graham 
[ 10 ]). 

The scaling law for time can also be derived from the scaled to actual 
distance ratio by noting that the time associated with the blast is the integration of 
the speed of the blast wave over the distance traveled. Thus, by Kinney and 
Graham [10], the actual time is related to the scaled time by 

. , . s (scaled _time) 

(actual _ time) = - (o.lo) 

ft 


and ft is the transmission factor for time taking also into account the ratio 
between the shock front speed, a, and the sound speed in the medium, ao. The 
expression for f is given as: 


f=fti — 
ci„ 


' Jf 

\Po J 


1/3 






1/2 


f n ^ 


yPoj 


1/3 




Kfj 


1/6 


(3.14) 


C. DYNAMIC LOADING OF BLAST WAVE ON AIRCRAFT 

1. Translational and Rotational Responses 

An object placed in the field of an explosive blast wave will experience a 
dynamic load characterized by the sudden increase in its peak pressure value 
followed by a gradual decaying process. The net effect depends on the 
distribution of the pressure wave across the aircraft’s surface, which depends on 
the orientation, geometrical shape and construction of the object, and its resistive 
forces such as bending moments, drag and aircraft’s thrust apart from its inertia. 
The end result is most likely a non-symmetrical distribution of the pressure 
intensity across the exposed area and this loading also changes with time. Thus, 
it is generally a difficult problem to solve analytically. 
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In this thesis, the main interest is to examine the translational and 
rotational responses of the aircraft as a whole. To accomplish this goal we 
assume the aircraft to be a rigid body with center of mass located at the origin of 
the aircraft coordinates (as in Chapter II except that we will expand the case into 
the third dimension). Consider the geometry of the three dimensional (3D) 
problem as shown in Figure 10. We assume the blast wave originates from a 
point source and diverges outward to the target. Different parts of the aircraft 
surface area exposed to the blast wave will then experience different peak 
overpressure intensities of different durations and different arrival times (because 
of their relative distances to the blast point). 
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Figure 10 3D geometry of the interaction of blast wave on the aircraft surface. 

A directed area element on the aircraft surface at point F can be described 
by the usual area element (A4) and normal vector, ni, describing the orientation 
of the aircraft surface. The point P is denoted by the vector a,-. A spherical 
diverging blast wave will travel in the radial direction along the directed line 
linking the blast point to the point where the pressure interacts with the target 
(point P). Denote this vector by Rbi, which can be derived if we know ai, and the 
distance between the center of mass of the aircraft and the blast point, denoted 
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by the vector, Rc. This distance can be taken to be the miss distance of the 
target (if we assume that the proximity fuse of the explosive munitions triggers at 
the point where the largest surface area of the target is sensed). The blast wave 
direction then follows along the direction of Rbi indicated by ui as 


Rsi _ 

\R,,\ \R^+a,\ 


(3.15) 


The elemental force acting at point P is then simply 


fi = Rc +«; l)A^(«,- -W/K- (3.16) 


where pi represents the instantaneous overpressure intensity of the blast wave 
given in (3.1). Note that pi is a function of time and distance between point P and 
the center of the explosion. 

It follows that — assuming the aircraft to be a rigid body — the total force, 
Fb, and moment, Mb, exerted by the blast wave can be expressed in terms of the 
sum of all the elemental force components over the exposed area of the aircraft 
body to the blast pressure wave: 

Fb= Z fi= Z Pi{t^R^+a^\)^A{n^■u^)u^ (3.17) 

Exposed _ area Exposed _ area 

Mb= Z Z K X A «/)«/] (3.18) 

Exposed _ area Exposed _ area 


Since momentum must be conserved in the process, the net effect of the 
application of these forces over a time increment t gives rise to the changes of 
both linear and angular momentum (acting on the center and about the center of 
mass) as follow 


F„ - Linear Resistance = — 

^ dt 

I (Fg - Linear_Resistance)dt = 


(3.19) 
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M, 


AngularResistance = Iq 


do) 

dt 


I (Mg - Angular_Resistance)dt = ] 


(3.20) 


where Linear Resistance and Angular Resistance represent the sum of all the 
elemental resistive force components. Putting (3.17) and (3.18) into (3.19) and 
(3.20) respectively and re-arranging them, we can express the variation in linear 
and angular velocity over time as 


(} Exposed _ area f 


- — f ( LinearResistance )dt 


(3.21) 


0) 


{t) = 0}^+^ ^ ^{a^x p.{t^R^+a^\)^A{n^■u.)u^\dt 

^ (} Exposed _ area f 


- —I (Angular_Resistance)dt 

t 


(3.22) 


The presence of the linear and angular resistive force components add an 
additional level of complexity when evaluating Equations (3.21) and (3.22). In 
general, these components vary with time as a result of the aircraft motion. 
However, even if the resistive component effects can be ignored, the analytical 
solution to (3.21) and (3.22) is only possible for simple geometrical aircraft 
shapes (in that we also have to assume that the orientation of the projected area 
does not change significantly with time due to the aircraft motion). Numerical 
solution (with the aid of computer simulation) may help to obtain a reasonable 
estimate of the linear and angular velocity variations for more general shape and 
different variations. Such simulations can also account for the resistive force 
components. 
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2. Approach to Numerical Simulation for ISAR Imaging 

ISAR imaging of an aircraft works by mapping the changes in the phase of 
the returned signals to the range and cross-range positions of the target 
scattering centers. Hence, our main interest is to derive from (3.21) and (3.22) 
the target’s radial range and angular position variations that are related to the 
phases of the returned ISAR signal for image generation. 

Consider the following 3D geometry of an ISAR looking at the aircraft as it 
is being hit by the blast wave. 



Figure 11 3D geometry of the ISAR staring at aircraft. 

From the geometry, the radial velocity detected by the radar can be given by 

v^(0 = v(0--^ (3.23) 

Hence, the range, R(t), and angular position, 6(t) of the target over the entire 
ISAR coherent integration period, T, can be obtained by integrating (3.23) and 
(3.22) over rto give 
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i?(0 = i?„ + Jv^a¥^ = i?o + J (3.24) 

e{t) = e^ + \\0){t)\dt (3.25) 

T 

From (3.24) and (3.25), we can obtain the range variation of a scattering 
point on the aircraft {Rp) located at P(x,y,z) (and represented by the vector P) by 
taking the magnitude of R(t) plus the projection of the vector P onto the fixed 
radar frame of reference (P is in the aircraft frame of reference and its frame is 
rotating with time varying 0 from the radar frame of reference). 

It is not trivial to provide a simple analytical simulation model for (3.24) 
and (3.25). However, if the distribution of R and 6 can be obtained over the 
integration period T and at intervals equal to the sampling time of the ISAR, S'„, 
through computer-based simulation, the results of this simulation can be used as 
inputs to the ISAR radar simulation model to give the solution for Rp (through 
coordinate transformation and projection). Following the previous discussion, Rp 
can then be used to map out the ISAR image through Fourier processing. Figure 
12 outlines this methodological approach. 




Figure 12 Approach to generate of target motion for ISAR simulation. 
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In the figure, the block on the left represents the simulation model for 
solving the equations of motion and provides a stream of R and 6 data over the 
period Twith interval to that of 5'„. The input to this computer model includes the 
target parameters (such as its location with respect to the radar, its shape and 
size, surface orientation, weight and motion), the blast wave parameters, its 
distance and direction from the target center of mass, the sampling time and the 
integration period from the ISAR. The output, in term of R and 6 together with the 
target scattering points’ locations, is fed into the ISAR simulation model on the 
right where they are used to generate the ISAR image. A distinct advantage of 
using this approach is that changes can be made to any one of the models with 
minimum influence on the others. 

D. BLAST WAVE EFFECT ON AIRCRAFT AND ISAR IMAGING 

In the development of Chapter II we showed that many ISAR systems rely 
heavily on the linearity of the target motions (in radial and angular velocity) over 
the integration period Tfor radar image generation. This is not desirable for use 
of ISAR in airborne target Battle Damage Assessment (BDA) since the target is 
likely to possess high dynamic motions as a result of the blast force interaction 
with its surface body. 

In Equation (3.21) and (3.22), we have shown that these motions contain 
higher order acceleration components that will alter the phases of the 
backscattered radar signal in a non-linear manner as it is collected by the radar 
receiver. Owning to the finite time duration of the blast pressure force, these 
non-linear components only exist as impulses for a fraction of the time equivalent 
to blast duration but their magnitude can cause significant changes to the aircraft 
motion. Also, these impulse forces leave behind residual velocity components 
for the remaining time of the radar coherent processing period, and are not easily 
mitigated by motion compensation algorithms. Hence, both the presence of non¬ 
linear and residual velocity components will act to interfere with the image 
processing algorithm and distort the ISAR image. 
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IV. TIME-FREQUENCY TRANSFORM METHODS OF ISAR 

IMAGING 

A. TIME-FREQUENCY REPRESENTATIONS OF SIGNALS 

The Fourier transform is a powerful and widely used method in signal 
processing and analysis because of its ability to decompose any arbitrary signal 
into a set of sinusoidal functions with different frequencies. The distribution of 
frequency components is useful for characterizing the properties and behavior of 
the signal, and on which signal manipulations such as filtering can be performed 
to alter them for specific application purposes. However, the Fourier operation 
and its inverse involve a one-to-one transitional relationship between the time 
and frequency domains. Consequently, the “information about the time 
localization of the frequency components” cannot be easily interpreted in the 
Fourier spectral domain^ (Hlawatsch and Boudreaux [6]). This limitation is not a 
problem for stationary types of signals but is generally true for time-varying 
signal. 

Many signals encountered in the real world have frequency content that 
varies with time. One common example is music where the harmonic content 
changes for different notes. Other common examples include biomedical 
signals, speech, vibrations and linear chirp pulses. The time-varying nature of 
such signal’s frequencies makes them rather difficult to represent in the time or 
frequency domain alone. The time-frequency Transform (TFT) helps to alleviate 
the situation by mapping a one-dimensional signal as a function of time into a 
two-dimensional time-frequency plane. Unlike the Fourier transform (which is 
always linear in nature), the time-frequency representation of the signal can be 
linear, bi-linear (quadratic) or can contain higher orders of non-linear terms 
depending on the method of signal representation. In most applications, the first 
two cases are widely employed and we shall discuss them as follow. 


^ Strictly speaking, this information is embedded within the phase of the complex valued 
frequency components but it is often awkward to operate on in the Fourier domain. 
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B. SHORT TIME-FREQUENCY TRANSFORM (STFT) 

1. Definition 

The short time-frequency transform (STFT) and its associated Gabor 
transform are the earliest forms of linear TFT’s and date back to the 1940s. In 
this case, the STFT can be defined as 

STFTitJ)=\sit')w*it'-t)e dt' (4.1) 

where s(t) is the signal and w*(t) is the complex conjugate of an analysis window 
function with a finite time duration. 

From (4.1), the STFT can be seen as taking the Fourier transform of the 
signal after it is multiplied by the complex conjugate of a shifted window function 
centered around t. Because the time width of the window function is finite and in 
general shorter than the signal duration, the effect is the “suppression of the 
signal outside a neighborhood around the analysis time point t’ = t." (Fllawatsch 
and Boudreaux [6]) The STFT represents simply the local spectrum of s(t’) at the 
analysis time t. 

The magnitude of the STFT(t, j) is often called the spectrogram of the 
signal and shows how the frequency spectrum varies as a function of time in the 
horizontal axis. This magnitude is often represented by a surface above the time- 
frequency plane and displays the spectral components at that time. 

2. Properties of STFT 

According to Fllawatsch and Boudreaux [6], the STFT also has a dual 
relationship in the frequency domain and can be expressed as 

STFTit, f) = j S{f W*{ff)e'"’'^''df ' (4.2) 

where S(f) and W(f) are the Fourier transform of s(t) and w(t), respectively. This 
relationship gives additional flexibility in representing the STFT of the same 
signal from the spectral domain. 
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There are other important properties of the STFT. The most significant 
property is that it satisfies the linearity principle where the STFT of the 
superposition of two signals can be represented by an equivalent linear 
combination of their respective signals’ STFT. 

^STFT{t,f) = c,STFT^{t,f) + c^STFT^{t,f) 

The linear behavior of STFT allows the decomposition of a complex signal into a 
basis set of any arbitrary functions to which the STFT can be easily applied. 

Other properties include the preservation of time/frequency under 
frequency/time shift: 

s{t) = ^ STFT(t,f) = STFT,(t, f - f,) (4.4) 

s(t)=s,(t-t,) => STFT(t,f) = STFT,{t-t,J)e-j2nft, (4.5) 


3. Time and Frequency Resolution 

The time and frequency resolution of the STFT are influenced by the 
choice of window function and, especially, its width. From (4.1), good time 
resolution can be achieved using shorter window duration. On the other hand, 
the dual relationship of Equation (4.2) tells us that a narrow bandwidth is required 
for good frequency resolution. But the uncertainty principle of time and 
frequency in the two domains does not allow the existence of an arbitrarily short 
time duration and small bandwidth window for the same signal (Hlawatsch and 
Boudreaux [6]). Flence, this principle inherently limits the time-frequency 
resolution of the STFT. 

Consider the following two cases (Fllawatsch and Boudreaux [6]). In the 
first case, we choose the window to be a Dirac delta function with perfect time 
resolution. Then, for the signal s(t), its STFT is given by 

STFT{t, /) = J s{t ')5{t (4.6) 
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The magnitude is essentially the same as the magnitude of \s(t)\ without any 
information about the frequency content of the signal. Similarly, if we let W(f) = 
5(f) and apply it to (4.2), the STFT of S(f) becomes 

STFTitJ) = j5(/= S{f) (4.7) 

which is the Fourier transform of s(t) - its magnitude bears no information about 
the time variation of the frequency components. 

The effect of time and frequency resolution can also be seen by referring 
to the illustrative example as follows. Figure 14(a) shows a signal whose 
frequency changes with time and Figure 14(b) is the Fourier transform of this 
signal. From Figure 14(b), we can see that the frequency domain does not 
provide any indication on how the signal frequency changes with time. 



(a) Time Domain Signal 



(b) Frequency Domain Representation 


Figure 13 Example of a time-varying signal with changing frequency. 
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The time-frequency representation of the same signal is shown in Figure 
14. This figure was generated using a sliding Gaussian window function. In 
Case 1, we set the width of the window function such that good time resolution 
can be obtained: notice the significant overlap in frequency. In Case 2, we 
broaden the window width to give good frequency resolution but in turn the time 
resolution suffers. 
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(a) Case 1: Small window size for good time resolution. 
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(b) Case 2: Large window size for good frequency resolution. 


Figure 14 Effect of Time and frequency resolution. 
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4. Shape of the Window Function 

The specific shape of the window function is important in reducing the 
side-lobe interference in the STFT representation. Typically, one requires the 
window function to taper toward zero smoothly and applies the Hamming, 
Hanning, Kiaser-Bessel, or Gaussian window function (Chen [2]). The best 
results have been shown to be achieved by Gaussian windows where the time- 
bandwidth resolution product is AtAf=l/2 (Refer to Chen [2]). The use of 
Gaussian windows in an STFT is known as a Gabor transform. 


C. BI-LINEAR TIME-FREQUENCY TRANSFORM 

Although the STFT form of signal representations has a nice linear 
property, the time-frequency resolution is limited by the uncertainty principle, 
which at best gives a time-bandwidth product of 1/2. Accordingly, Cohen [12], 
Hlawatsch and Boudreaux [6] and many others have shown that significant 
improvement in the time-frequency resolution can be obtained by representing 
the signal in the quadratic or bi-linear TFT form. The quadratic representation is 
also “an intuitively reasonable assumption when we want to interpret” the time- 
frequency of the signal as a form of energy distribution (Hlawatsch and 
Boudreaux [6]). 

To classify a TFT signal representation in the bi-linear form, it must satisfy 
the following marginal properties related to the energy density. 


\TFT{tJ)df = Pit) = \s{t)\" 
\TFTit,f)dt = Pif) = \Sif){ 


(4.8) 


where p(t) and P(f) are the instantaneous power and the spectral energy density 
of the signal s(t) respectively. As a consequence, the signal energy E = j \x(t)fdt 
= 1 \X(t)fdf can be derived by integrating TFT(t,f) over the entire time-frequency 
plane. 
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The most difficult part in interpreting the bi-linear TFT representations is 
the presence of the cross-interference terms which occur as a result of the 
quadratic superposition principle. Consider the spectrogram of the STFT (which 
is often loosely interpreted as a bi-linear TFT representation although it does not 
meet the marginal properties according to Fllawatsch and Boudreaux [6]). 

SPEC{t, f) = \STFT{t, /)|" (4.9) 

We can easily show that, for the signal s(t) = si(t) + S 2 (t), the spectrogram of this 
composite signal is not SPECi(t,j) + SPEC 2 (t,j) but contains cross interference 
terms: 


SPEC{t, f) = \STET^ (t, f) + STET^ (t, f)f 

= \STET^(t,f)f + \STFT^(t,f)f + 2STFT^(t,f)STFT^(t,f) (4.10) 
= SPEC^ (t, f) + SPEC 2 (T /) + Cross _ Intereference 


This behavior is generally true for all types of bilinear TFT representations 
although the magnitude of the cross term varies with different types of bi-linear 
TFT transforms. Also, the more complex the signal is, the more the linear 
combination of its basis set is used to represent the signal; hence, the more 
cross interference terms that will appear in the spectrogram and bi-linear TFT 
representations (Fllawatsch and Boudreaux [6]). 


There are many forms of bilinear TFT. The most basic is the Wigner-Ville 
Distribution (WVD), which is defined as the Fourier transform of the time- 
dependent autocorrelation of the signal R(t, t’) (Chen [2]) 


WVD{t, /) = j R{t, t ' 


where R(t, t’) is chosen as 


R(t,t') = s 



( 



s 

\ t - 

5* 

t - 


1 2; 


1 2j 


(4.11) 


(4.12) 


Because of this relationship with the auto-correlation function, the WVD is 
known to have very good time-frequency localization and hence good time and 
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frequency resolution. But this distribution is also known to possess large cross 
interference terms that hinder its direct usefulness for signal analysis and 
processing. 

A more generalized form known as the “Cohen class” of transforms 
extends the use of the WVD to include other members of the bilinear TFT 
representation with a kernel function (l)(t, t’). The general form of the Cohen class 
is given by (Chen [2]) 



(j){t-uX)e dudt' 


(4.13) 


Note that the WVD is simply a subset of the Cohen’s class with t’) = 5(t). The 
significance of the Cohen class is that different types of kernel function can be 
designed to reduce the cross-term interference of the bilinear TFT representation 
but, in the process, compromise the time-frequency resolution (compared to 
WVD — which is known (Chen [2]) to possess the best time-frequency 
resolution). Two well-known distributions in this category are the Choi-Williams 
distribution (CWD) and the cone-shaped distribution (CSD). Further 
development of these topics is beyond the scope of this thesis, and we shall not 
discuss them further. For detailed reviews, refer to excellent texts by Chen [2], 
Fllawatsch and Boudreaux [6] and Cohen [12]. 


D. TIME-FREQUENCY BASED IMAGE FORMATION 

In Chapter II, we noted the stringent radial and angular motion 
requirements imposed on ISAR systems for high quality image generation: both 
velocities must essentially remain constant throughout the coherent integration 
period. The constant radial velocity allows the radar to accurately estimate the 
target’s linear motion and compensate for it accordingly. To maintain a simple 
direct relationship for mapping the returned phase variations to the cross-range 
position of the target’s scattering centers, the angular rotation rate must be linear 
and within an upper bound such that the viewing angle of the ISAR is small 
enough to prevent cell migration. 
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However, when the aircraft is being hit by explosive munitions, the blast 
force will cause changes in both linear and angular momentums. The net effect 
is that the target will exhibit complex variations in its radial and angular velocities, 
occurring over a very short time interval (typically, msec) and well within the 
coherent integration time (usually measured in sec). Most translational and 
rotational motion compensation techniques are inadequate to provide good range 
and Doppler tracking of these complex motions. As a result, the image becomes 
distorted when conventional Fourier-based techniques are employed for image 
construction. 

The complex variation of the radial and angular velocities from the blast hit 
can be viewed as a time-varying frequency shift of the returned signals within the 
integration period. Although these signals pose problems for conventional FFT 
processing in ISAR systems, the signals do contain information about the 
dynamic behavior of the aircraft as it is hit, and this information can be extracted 
using TFT processing. 

From an ISAR system architecture point of view, the front end processing 
remains largely the same. The SFW returned signals are detected and sampled 
as I- and Q- data, and shipped off to the appropriate range and Doppler tracking 
algorithm for motion compensation. Synthetic range profile processing is then 
performed on Vc\e Mx N complex-valued sampled data array to obtain the range 
information as usual. The result is also diM xN complex-valued matrix array, H. 
Note that, in accordance with our earlier definition, each row of the matrix 
(denoted as h( m , 1:N ) where m = 1, 2, M ) represents a time-history series 
comprised of N pulses of range cells. Each column of the matrix (denoted as h( 
1:M, n ) where n = 1, 2, N) represents a collection of time-histories of the 
target from M burst at the particular range cell index n. 

In conventional ISAR image processing, the direct Fourier transform 
(usually implemented by an FFT) is applied to each column h( 1:M, n) \.o obtain 
the cross-range position of the target scattering points in that range cell index. 
The Fourier processing still yields a one dimensional column vector: 
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-.M ,n) = FFT„ {hil'.M ,n)} (4.14) 

This process is repeated for all N range cell indices. The result is still a 2D MxN 
matrix array with the magnitude of each element indicating the presence of the 
target point scatterer in the particular range and cross-range positions. 

TFT’s differ from the conventional Fourier transform method of processing 
in that the windowed transform (can be linear STFT or bilinear TFT) maps the 
one-dimensional collection of time-histories of M bursts in each column h(:, n ) 
into a two-dimensional M x M plane revealing information about the Doppler 
variation along the time-histories in that cell index. 

,\:M ,n) = TFT^{h{\-.M ,n)] (4.15) 

This process is repeated for all range cells and we get a three dimensional 
complex-valued array of MxMx A^Doppler-time-range matrix. Figure 15 depicts 
the image processing based on TFT method. 
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Figure 15 ISAR image processing based of TFT method. 


A series of M ISAR images can then be obtained from the three 
dimensional matrix array by summing over all the range cell indices in each m 
time-history point as shown in the figure below. 
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Figure 16 Extraction of ISAR Image from the Doppler-time-Range Matrix. 


The image extraction process can be summarized by the following expression: 

I^{l:M,l:N) = D{l:M,m,l:N) for m = 1, 2, M (4.16) 

When the same process is repeated for each time-history series, the time-history 
series display of the 2D ISAR image with indices running from m = 1 to M is 
analogous to a moving picture of the aircraft within the integration period (as 
opposed to the static one by the conventional Fourier-based construction 
method). Flence the dynamic motion of the aircraft due to the blast effect is 
revealed by these sets of ISAR images. 

The image resolution from TFT’s is generally coarser in comparison to the 
image obtained by conventional FFT. The resolution also depends on the width 
and shape of the window function for linear TFT. Bilinear TFT’s will give better 
resolution, but the presence of cross interference terms may act to distort the 
image. 
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V. DEMONSTRATION OF BLAST EFFECT 
CAPTURINGBY ISAR USING TIME-FREQUENCY 

TRANSFORM 

A. BLAST EFFECTS AND THE ISAR SIMULATION MODEL 

To demonstrate the possibility of using Time-Frequency Transform (TFT) 
methods to ‘see’ the dynamic behavior of an aircraft target in ISAR imaging, a 
simple simulation model was developed as part of this thesis. This simulation 
contains two primary sections following the methodological approach described 
in Chapter III. The first section aims to simulate the response of the aircraft as 
the result of a blast hit. The second part of the model extends the usual ISAR 
imaging simulation models to include a time-frequency method of radar imaging 
based on the work of Chen [2]. 

The simple model developed for this demonstration does not provide a 
complete description of the real world problem at hand involving target BDA. 
There are several limitations associated with this model that deserve comment: 

1) Only two-dimensional geometries are considered to ease the 
development and thesis illustration. Nevertheless, one can extend the 
model to the three-dimensional case with careful considerations about the 
target roll, pitch and yaw effects, and with some manipulations of the 
frame of references using matrix transformation and vector operations. 

2) Initially, we shall consider only the translational response of the air 
target in this demonstration. The model’s implementation requires 
analytical expressions for translational and rotational velocities. It is 
easier to obtain the translational motions from (3.21) with some 
reasonable assumptions (as we shall see later). Flowever, the angular 
dependence in Equation (3.22) contains terms related to the geometrical 

6 Some of these models are easily available in or as part of publications related to ISAR 
imaging in the courtesy of the author (s). For this thesis, we would like to acknowledge that our 
model is a modified version of the Matlab ISAR simulation program using conventional FFT 
processing from Jae, Thomas and Flores [3] to help in hastening the model development time 
compared to starting from scratch. 
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distribution of the projected aircraft area and the distribution of the 
pressure wave along its surface that needs to be properly accounted for in 
order to obtain reasonable estimates of the net effect in angular velocities. 
There is no simple analytical solution to this problem, and it requires a 
large scale numerical approach which is beyond the scope of this thesis. 

3) Linear time-frequency transform methods (in particular the Short 
Time-Frequency Transform or STFT) will be used in this study. The main 
concern with linear TFT’s is the time-frequency resolution, which (at best) 
gives a time-bandwidth product of only y 2 . We will have to evaluate its 
adequacy for our purpose before we proposed further improvements to it. 

The code for the simulation program was developed using Matlab 6.5 with 
time-frequency toolbox provided for by Auger et al [7] & [8] available at 
http://crttsn.univ-nantes.fr/~auqer/tftb.html , Aug 2004. The entire code is given in 
Appendix A. 

B. SIMULATION OF AIRCRAFT RESPONSE TO BLAST WAVE 
INTERACTION 

The primary goal of this part of the simulation is to obtain the translational 
motion response of the aircraft from the blast hit. This response appears as a 
form of variation in range as a function of time and can be derived from the 
equation of motion for translational motion (3.21) and the blast wave expressions 
given in (3.1). There are certain assumptions needed in order to simplify the 
derivation. They are described as follows. 

1. Blast Pressure Distribution on Aircraft 

To ease the stringent geometrical requirement of the non-symmetrical 
aircraft surface, we assume the blast wave to be a plane wave with equal 
pressure distribution on the exposed area of the aircraft body. The exposed area 
is also assumed to remain constant with time during the evaluation period. By 
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ignoring the resistive force components, Equation (3.21) can be re-arranged and 
simplified to become 


v{t) = v^+^{^p{t,R^)dt Yj ^ 


M, 


G \ l 


{n. ■ Ui)Ui 


= V +- 


M. 


G 


J p(t,RJdt 

K t ) 


J y Exposed _ area J 

U: 


(5.1) 


where Vo is the initial velocity, Mg is the mass of the aircraft, Ap is the projected 
area of the aircraft exposed to the blast pressure given by 




\^Exposed _area J 




(5.2) 


and p(t, Rc) is the instantaneous pressure wave described by Equation (3.1). Rc is 
factored into the original p(t) since the parameters determining the shape of p(t) 
— Po, ta, td and a — are actually dependent on the blast distance, Rc. The unit 
vector Hi represents the direction of the net blast effect which is normal to the 
aircraft surface subjected to the blast. 


Now that the blast pressure acting on the aircraft body is independent of 
its direction and that the distance Rc can be taken to have constant value under 
the planar wave assumption, the impulse response (which is integral of p(t) over 
time t) can be evaluated separately. If we denote the impulse response of the 
target by I(t), then 


l{t) = ^ p{t\R^)dt' 
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tARc) 
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dt ' where 


where po(Rc), ta(Rc), td(Rc) and a(Rc) can be expressed as a function of the blast 
distance using the empirical Equations (3.2), (3.4), (3.6) and (3.8) respectively. 

We denote po’ as po(Rc), ta’ as ta(Rc), td’ as td(Rc) and a’ as a(Rc). Applying 
the appropriate boundary conditions, we can obtain the solution for I(t) as 
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2. Geometry of the Problem and its Analytical Solution 

We are interested in the radial velocity component since this is detectable 
by the radar as a Doppler shift. To obtain the radial velocity variation, consider 
an exaggerated two-dimensional geometry of an ISAR looking at the aircraft 
while an explosive blast occurs near it as shown in Figure 17. 



Figure 17 2D Geometry between ISAR, aircraft and blast point. 

For the geometry of the figure, the force acting along the radar line-of- 
sight (LOS) is given by 

(5-5) 
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where Ap is the projected surface area of the aircraft to the blast, p(t, Rc) is the 
instantaneous blast pressure wave as in (3.1), Rc is the blast distance, and 6b is 
the direction of the blast wave. 

We use the following sign convention: the radial velocity is positive when 
the target is approaching the radar and negative otherwise. Hence, the resultant 
radial velocity variation can be simply obtained using conservation of momentum 
(ignoring all other resistive forces) as follows 

-J Flos [^r (0 “ v,(0)] (5.6) 


Re-arranging (5.6) and substituting (5.5) into it, we get 

A„ sin 


vM) = vM- 
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A^ sin 
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Mr- 


-j p(t,RJdt 


(5.7) 


(Equation (5.7) can also be obtained from the projection of (5.1) onto the radar 
LOS and the integral is exactly the same as in Equation (5.4).) 

For the ISAR simulation model, the target linear motion is expressed as 
the range variation R(t). To obtain R(t), we simply integrate v/t) with respect to 
time using the following kinematic relationships. 

R t 

R{t) = R{Q) + \v^{t')dt' (5.8) 

t 

^„sin6*„ f 

= i?(0) + v^(0)r- " ^ \mdt' 

t 


The result of the integral l(t) is expressed below, which represents the 
translational dynamic behavior of the aircraft due to blast effect. 
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3. Aircraft Profile 

Following Equation (5.2), we can define the overall projected area of the 
aircraft exposed to the blast as the summation of the projected frontal area of its 
body, Af, and the projected side area of its body. As, and expressed in term of the 
aircraft’s orientation with respect to the blast direction as in Figure 17. 

A^ = \Ap cos{0g +0)\ + \Ag sin(^^ +0)\ (5-11) 

The front and side profiles of the aircraft used in this simulation have a 
shape shown in Figure 18. Flence, Ap and As are given by 

AF=7r — + {W^-h)h^ ( 5 . 12 ) 

Ag = hL + Ljhj 
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Figure 18 Aircraft area profile. 

4. Blast Parameter 

Peak overpressure po, arrival time ta, blast duration td and the waveform 
parameter a are needed to fully describe the blast wave as a function of time, po, 
td and a can be easily obtained using the empirical equations given in Chapter III 
— (3.2), (3.6) and (3.8) respectively. The arrival time 4, however, requires 
integration of the speed from the center of explosion to the distance the shock 
front hits the aircraft. To simplify the simulation, tabulated forms of these data 
are extracted from Kinney and Graham [10] where they are scaled to 1kg TNT 
equivalent values at standard atmospheric pressure and temperature conditions. 
The actual values can be obtained by using Hopkinson’s scaling law described in 
Chapter III. 

For our simulation, the values for the explosive parameters are: explosive 
weight = 35kg with a TNT equivalent ratio of 1.5 for FIE type of explosive. These 
are figures appropriate to many long range Surface-to-air and Air-to-air missiles. 
The blast distance is chosen to be 5m, which is within the miss distance of most 
missiles today. 
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C. SIMULATION OF ISAR IMAGING USING TIME-FREQUENCY 

TRANSFORM 

The simulation model for generation of the ISAR image consists of three 
distinct parts: (1) Target model representation, (2) Radar signal generation and 
(3) ISAR Image processor. The description of each of these portions is as 
follows. 

1. Target Model Representation 

It is a common practice in most ISAR simulations to represent the ISAR 
target as a series of point scatterers outlining the shape of the target. The 
number of scattering points and the reflectivity of each point can be arbitrary 
chosen as long as they are within the image plane. For our case, 49 scattering 
points were used to outline a target and their reflectivity set to one (i.e., uniform 
signal reflection). 

In the model, the target is free to assume any initial orientation with 
respect to the radar by varying Oq. The motion for each of the target scattering 
points is modeled in term of range variation using Equation (2.4), and since the 
target is treated as a rigid body, the range variation for all the points is influenced 
by the same kinematic equations for R(t) and 6(t). The expression used for R(t) 
follows from the discussion of section (5.8), and 6(t) is assumed to behave as in 
Equation (2.3). 

2. Radar Signal Generation 

The image plane for the ISAR processing is represented by 64 x 64 cell 
indices — a N x M matrix array (note that the range and cross-range axis are 
swapped in this presentation). The radar parameters,/c, Af, N pulses, M burst, 
pulse interval 1/ PRF and integration period T together with the initial range Ro 
and rotation rate coo are specifically chosen to optimize the overall radar system 
such that an ideal image can be obtained without the need to consider range 
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alignment, range and Doppler tracking, and motion compensation in the process. 
A summary of these parameters are given in Table 1 below. 

The radar signal generation follows the theory in Chapter II and is simply: 
compute the radial range R(t), angular variation 6(t) and stepped frequency 
increases/„,„ under (2.8) at each sampling time tm.n in (2.9) for T, and re-arrange 
them as a N x M matrix. These three time-varying matrices are imposed as 
phases onto the target model and the signal is taken as the sum of these phases 
from all the scattering points according to (2.10). 


Table 1. Radar Parameters 


Parameters 

Values 

Remarks 

Target Initial Rotation Rate, o 

3.77deg/s 

Optimized for the target length. 

Target Initial Range, Rq 

22.567km 

Optimized for range alignment. 

Carrier Frequency,yc 

3GHz 

Typical Operating value. 

Stepped Frequency, Af 

2.2MHz 

Optimized for range resolution. 

Number of Pulses, N 

64 

Arbitrary chosen. 

Number of Burst, M 

64 

Arbitrary chosen. 

Pulse Interval, T 2 

150|isec 

Optimized for Doppler resolution. 

Pulse Width, Ti 

0.454|isec 

Ti=l/Af 

Integration Period, T 

0.6144sec 

T = NxMxT2 

Max Target Length, x 

30m 

x=2(Arcf/A 

Range Resolution, Ar 

1.06m 

Ar=c/(2NAf) 

Cross-Range Resolution, Arc 

1.23m 

Arc=X/(2a)T) 


65 






3. ISAR Image Processor 

The image processing of the radar signal begins by generating the 
synthetic range profile of the target according to Equation (2.11). Subsequently, 
the resultant data is passed to another FFT process to construct the conventional 
two-dimensional mapping of the target image using (2.12). 

On the other hand, time-frequency image construction follows the 
procedures discussed in Chapter IV: the STFT is applied to each row of the Nx 
M synthetic range profile matrix as in Equation (4.15). The result is a three- 
dimensional Doppler-time-range mapping of the target. Each of the M time- 
histories is then extracted for display accordingly as shown in Figure 16. The 
operation is carried out using Equation (4.16). 
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VI. RESULTS AND DISCUSSIONS 


A. TIME-FREQUENCY REPRESENTATION OF ISAR IMAGE 
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(a) Target model. (b) ISAR image. 

Figure 19 Target model representation and its ISAR Image (from FFT process). 

Figure 19(a) above show the model of an aircraft used in this study. The 
model is represented by the series of 49 point radar scatterers spread over the 
body length of 35m and wing span of 20m. Figure 19(b) is an ISAR image of the 
aircraft model generated using the conventional Fourier process (in particular 
Fast Fourier Transform, FFT) described in Chapter II. 

When the same target model is analyzed using the Short Time-Frequency 
Transform (STFT) method (discussed in Chapter IV), a successive series of 64 
time-history ISAR images can be displayed instead. Each of these images shows 
the same outline of the aircraft image. Figure 20 below shows an example of 16 
images extracted from the time-history series oi m = 33 to 48 (at time t between 
0.3168 to 0.4608sec where t = 0 indicates the beginning of the coherent processing 
period). 

The image resolution is comparatively poorer — because of the tradeoff in 

frequency resolution given to the time resolution. Conventional FFT’s operate 

over the entire time-history series of M cells (in our case M = 64) but produce 

67 




only a one-dimensional mapping in the spectral domain. Whereas in STFT, the 
Fourier processing is carried out with a sliding time window function to yield a 
two-dimensional time-frequency mapping of each time-history series m where m 
= 1, 2, M. The width of the window is generally chosen to be shorter than M to 
give better time-resolution, which is important (as we shall see later) for 
observing the dynamic behavior of the target within the processing period. (In its 
extreme, we can show that exact FFT image is obtainable from the STFT simply 
by setting the window width to cover the entire M cells). But the instantaneous 
frequency resolution suffers as a result, and this frequency is directly related to 
the cross-range resolution of the ISAR image. 
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Figure 20 Sample of ISAR from STFT. 
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B. TRANSLATIONAL RESPONSE OF AIRCRAFT DUE TO BLAST 

EFFECT 

1. Blast Characteristic 

In this study, we place the blast directly between the target and radar 
along its line-of-sight (LOS) so that the magnitude of translational response will 
be maximal in the radial direction. In the first case studied, we only let the side 
of the aircraft be exposed to the blast wave. The aircraft has a side area of 64m^ 
and weight of 22tons. The translational responses in term of acceleration, 
velocity and range variation are as shown in Figure 21 and Figure 22. 



Figure 21 Target acceleration and velocity variation (Blown up version). 



Figure 22 Target range variation. 
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As expected, the overall acceleration directly follows the impulse response 
of the blast pressure wave. The presence of this acceleration force gives rise to 
the non-linear variation in the radial velocity which eventually dissipates away. By 
conservation of momentum, the change in velocity will cease as the force is 
taken off the load and is set to a new (constant) value. We note also that the 
range variation from zero to 3m results from the velocity components. Since 
these effects happen within one integration period, the ISAR system will not have 
enough time to compensate for them. Hence the blast effects are effectively 
introduced as errors in the imaging process. 

2. ISAR Image from Conventional FFT 

These errors are manifest as distortions to the radar image as shown in 
Figure 23(b) when constructed using a conventional FFT. Notice the ‘smearing’ 
of the scattering points due to the migration of the range and cross-range cells in 
the presence of the radial velocity component, and the creation of a seemingly 
‘double image.’ We shall discuss this later when we examine the STFT image. 
But other than the presence of cell migration, the ISAR image does not reveal 
anything about the target’s behavior within the integration time. 




10 20 30 40 50 60 

Cross-Range Bin 


(a) Undistorted. (b) Distortion due to blast effect. 

Figure 23 ISAR image. 
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3. ISAR Image from STFT 

The ISAR images constructed using STFT are shown below in Figure 24. 
The period of interest is around the time when the blast interaction occurs, which 
is the first O.Ssec. In the images, we can observe the sequence of events that 
happen to the aircraft when it is being hit by the blast wave. During the first 
O.lSsec, the position of aircraft image is the same as before. As the blast wave 
arrives, the aircraft image begins to blur because of the velocity variations. This 
happens for about the same duration as the blast wave interaction with the 
aircraft body. When the variation stops, the aircraft image becomes gradually 
clearer (at t = 0.2sec) but seems to be offset to a new position, which is an 
expected effect of a result of the uncompensated velocity components present in 
the image processing. 


TH1 (t = 0 0096sec) 
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TH9 (t = 0 0864sec) 
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TH16(t = 0.1536sec) 



(a) At time t = 0 to 0.1536sec. 
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(b) At time t = 0.1536 to 0.3072sec. 


Figure 24 STFT ISAR images for observing the blast effect. 

From these STFT images, we also make two interesting observations: 

1) There is a direct correlation with these images to the one obtained 
using FFT. We notice that the ‘smearing’ of the image earlier is actually 
caused by the ‘excessive’ cell migration from the nonlinear variation of the 
radial velocity when the blast wave hits the aircraft. Also, the brighter 
‘double image’ is simply the result of continuous imposing of the offset 
aircraft image after the blast. Flowever, note that the blast has unusually 
large effect along the cross-range axis rather than the range (The cross¬ 
range cell shift is counted to be about 8-10 cells whereas the range cell 
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variation is only 3 cells) although the only velocity component is presence 
in the radial direction. 

2) The image resolution is much poorer compared to Figure 20. This 
is a consequence of reducing the size of the window function to only y4 of 
the M cells so that the time resolution is appropriate to reveal the blast 
event happening around the time t = 0.165 to 0.19sec. However, we are 
unable to capture the exact time that the event occurred. Notice that the 
image blurring begins at t = 0.1344sec and stops only after t = 0.2016sec. 
This is caused by the overlapping of the edges of the window functions at 
the times when the window is centered in the period that the blast occurs. 
The time resolution can be improved with shorter window width or sharper 
cutoff at the edges but such improvement is at the cost of poorer 
frequency resolution. 

These two points shall be discussed in more detail in the later section. 

We verify the consistency of the result obtained with a second case in 
which the orientation of the aircraft is set to face away from the radar and only its 
back area is subjected to the blast wave. Consequently, the magnitude of the 
radial velocity variation is reduced from 5.7m/s to less than 0.7m/s. (But the blast 
characteristic remains the same.) 
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Figure 25 Acceleration and Velocity variation for 2"^ case. 
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Figure 26 ISAR image from FFT for 2^^ Case. 
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TH7 (t = 0 0672sec) 



TH11 (t = 0 1056sec) 



TH15 (t = 0 144sec) 



TH4 (t = 0 0384sec) 



TH8 (t = 0 0768sec) 




TH12 (t = 0.1152sec) 



TH16 (t = 0 1536sec) 



(a) At time t = Oto 0.1536sec. 


74 



















TH17 (t = 0 1632sec) TH18 (t = 0 1728sec) TH19 (t = 0 1824sec) TH20 (t = 0 192sec) 









TH21 (t = 0 2016sec) TH22 (t = 0 2112sec) TH23 (t = 0 2208sec) TH24 (t = 0 2304sec) 





TH25 (t = 0 24sec) TH26 (t = 0 2496sec) TH27 (t = 0 2592sec) TH28 (t = 0 268Bsec) 





TH29 (t = 0 2784sec) TH30 (t = 0 288sec) TH31 (t = 0 2976sec) TH32 (t = 0 3072sec) 





(b) At time t = 0.1536 to 0.3072sec. 


Figure 27 STFT ISAR images for 2'^^ case. 


Similar cell migration effects are observed in this case, except that the 
extent of the shift is slightly lessened compared to the first case. This behavior is 
expected since the velocity variation is lowered in the second case. The cell 
range-migration is hardly noticeable (virtually zero) whereas the cross-range cell 
variation is about 6-8 cells. Also, the image blurring is less significant in 
magnitude and duration because of the smaller blast force and shorter blast 
duration. Other cases were also simulated and examined with different target 
orientations. Their results, and the associated conclusions, do not differ from 
those described above. 
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C. SELECTION OF WINDOW FUNCTIONS 

In the above examples we have seen the important relationship between 
the time-domain and frequency-domain resolutions (a relationship which, 
unfortunately, is directly linked to the image quality). For STFT, these resolutions 
are determined by the size and shape of the window function, and the uncertainty 
principle associated with the Fourier transform. 


1. Resolution of Rectangular and Gaussian Window Functions 

In general, there are numerous window functions to choose from but, for 
illustration purposes, we will restrict our attention to the two cases for linear 
TFT’s: the simple rectangular function; and the Gaussian function. The 
rectangular function can be expressed as 


w{t) = 


, a a 

where - — — 

2 2 


|1 

0 otherwise 

W (/) = a sin c{naf) 


( 6 . 1 ) 


where a is the width of the window. And the Gaussian function expression is 
given by 


w(f) = 



W(/) = (2af^7r''V 2 ^ 


( 6 . 2 ) 


The resolutions for the two window functions — in both time and 
frequency — can be obtained using (Chen [2]) 


At = 


J {t- Mf\w{t)'\ dt 

-00 

00 

J \w{t)f dt 


,A/ = 


]{f-lii,f\W{f)fdf 

-00 

00 

\\w{fidf 


1/2 


(6.3) 


where 
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\t\w{t)\"dt \f\W{f)fdf 

Mt==l -- 

\\w{t)\"dt \\W{f)fdf 

-00 -00 

Thus, putting (6.1) into the Equations in (6.3), we derive the time and 
frequency resolutions for the rectangular function to be At = a, Af = 1/a, and its 
time-bandwidth product AtAf = 1. Both resolutions are directly related to the 
width of the window and multiplicative inverses of each other. 

Also, putting (6.2) into (6.3), we obtain for the Gaussian window function. 
At = a/v^, Af = 1 / (av^) and AtAf = Vi. The two resolutions are related to the 
shape of the window function, which is determined by cr. This gives us an added 
advantage of stretching the window width to cover the entire time-histories of M 
cells. 


2. Effect of Using Rectangular Window Function in STFT Imaging 

In our earlier demonstration (the first case shown in Figure 24), we used a 
rectangular window function with a window size a = V 4 of 64 cells = 16 cells. The 
interval between each cell was N x T2 = 64 x ISOjusec = 9.6msec. By setting the 
window width to 16 cells, we obtained the resolutions At = a = 153.6msec and Af= 
1/a = 6.51 Hz. 

For that example, the blast event happened between time steps 0.165 and 

O.lQsec and was only 25msec in duration. Thus, the selected time resolution is 

about six times larger than the blast duration. This fact explains why the STFT 

ISAR image could not faithfully capture the event. The reconstructed images 

began to mark the blurring effect at t = 0.096sec and end at about t = 0.2496sec. 

The difference of 152.6msec is equivalent to the time resolution of the rectangular 

function. To obtain accurate tracking of the event, the window size should be 

smaller than 25msec, which (for our case) is one to two m-cells. But in doing so, 

the frequency resolution will decline significantly to about 1/a = 1/ (2 x 9.6msec) = 
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52.1Hz. At this frequency resolution, it is virtually impossible to map out the image 
at all. Below is an example showing the serious degradation of the image by 
setting the window size to 2 m-cells. 


TH1 (t = 0 0096sec) 



TH13 (t = 0 1248sec) 



TH2 (t = 0 0192sec) 



TH10 (t = 0 096sec) 



TH14 (t = 0 1344sec) 



TH3 (t = 0 0288sec) 



TH11 (t = 0.1056sec) 



TH4 (t = 0 0384sec) 



Figure 28 STFT ISAR Time-history images with width set to 2 m-cell. 


The Doppler resolution for image quality is inversely proportional to the 

integration time T. Since T = 0.6144sec, we obtain the required frequency 

resolution to be 4 /d = i / T = 1.62Hz. But the frequency resolution of the above is 

four times coarser than Afo. As a result, the image quality is generally poorer 

compared to the conventional FFT ISAR imaging. We can improve the image 

quality by extending the window size to y 2 of the M cells (or more) but that would 

mean losing the time resolution and the dynamic behavior of the aircraft 

response could no longer be tracked. An example of this effect, with the window 

size set to 55 m-cells, is shown in Figure 29. Notice that the position offset is 
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shown right at the start of the integration time while the blast effect occur at t = 
0.165sec. But the quality of the images is better than previous ones. 


TH1 (t = 0 0096sec) 


TH5 (t = 0 048sec) 


TH2 (t = 0.0192sec) 



TH6 (t = 0 0576sec) 


TH3 (t = 0 0288sec) 


TH7 (t = 0 0672sec) 







TH9 (t = 0.0864sec) 



TH13 (t = 0 1248sec) 




TH4 (t = 0.0384sec) 



TH8 (t = 0 0768sec) 



TH12 (t = 0 1152sec) 



Figure 29 STFT ISAR Time-history images with width set to 55 m-cell. 


3. Effect of Using Gaussian Window Function in STFT Imaging 

Compared to the rectangular window function, the time-bandwidth product 
of the Gaussian window of only y 2 allows better time and frequency resolution to 
be obtained with lower limits on each. Figure 30 shows a graph comparing the 
Gaussian and rectangular window function in term of their time versus frequency 
resolution. The graph is obtained using their respective time-bandwidth product. 
From the graph, we can easily see that the Gaussian window requires a smaller 
window size to obtain quality image frequency resolution. Conversely, for the 
same time resolution it has better frequency resolution than the rectangular 
function. 
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Figure 30 Time and frequency resolution comparison chart for rectangular and 

Gaussian window functions. 


TH1 (t = 0 0096sec) 



TH5 (t = 0 048sec) 



TH9 (t = 0,0864sec) 



TH13 (t = 0 1248sec) 



TH2 (t = 0 0192sec) 



TH14 (t = 0 1344sec) 



TH3 (t = 0.0288sec) 



TH11 (t = 0,1056sec) 



THIS (t = 0 144sec) 


- V 


TH4 (t = 0.0384sec) 



TH12 (t = 0 1152sec) 



TH16 (t = 0 1536sec) 



(a) At time t = Oto 0.1536sec. 
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TH17 (t = 0 1632sec) 



THIS (t = 0 1728sec) 



TH19 (t = 0 1824sec) 



TH31 (t = 0 2976sec) 


■ft 


TH20 (t = 0 192sec) 



TH24 (t = 0 2304sec) 


TH28 (t = 0 2688sec) 



TH32 (t = 0 3072sec) 


S3> 


(b) At time t = 0.1536 to 0.3072sec. 


Figure 31 STFT ISAR Time-history images with Gaussian window function. 

Figure 31 shows the same 32 STFT images generated using the 
Gaussian window and should be contrasted with the rectangular window used 
earlier. To enable a direct comparison, we set the time resolution of the Gaussian 
window to be the same as in the first (rectangular window) case i.e. 16 x 9.6msec 
= 153.6msec. Since the time and frequency resolutions are determined by the 
shape of the window function rather than its width, we use a= Atv^ = 0.2172 to 
obtain the corresponding frequency resolution as Af = 1 / (av^) =3.26Hz. This 
value is only half of the previous case. 
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We can see that the time resolution remains the same by observing the 
blurring effect starting at t = 0.096sec and ended slightly early at t = 0.2208sec. 
The difference of about 124.8sec is close to the time resolution. However, the 
images’ quality has improved over that of the earlier case (as a result of the finer 
frequency resolution). 


D. EFFECT OF RADIAL VELOCITY VARIATION ON ISAR IMAGE 

In the images obtained using FFT or STFT, cell migration in both range 
and cross-range directions was quite noticeable — particularly along the cross¬ 
range axis. This migration caused the ‘smearing’ of the image and the blurring 
effect. 

To quantify the extent of the cell migration along the range due to the 
presence of the radial velocity component, we use the expression (2.27) given in 
Chapter II to obtain a plot of AN versus the radial velocity (Figure 32). For a 
change in the radial velocity that corresponds to about 4-5 m/s in our 
demonstrated blast characteristic, the estimated shift from Figure 32 below is AN 
^4 range cell. This is about the same amount of range cell shift counted in the 
ISAR images (which is 3 cells given in the previous section). Surprisingly, it is 
not a significant contribution to image distortion but only causes a slight 
‘smearing’ along the range cell seen earlier. 



Figure 32 Extent of cell migration due to velocity. 
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Comparatively, cell migration along the cross-range is much higher 
although there is no angular variation components present in the demonstration. 
To explain this strange effect, consider the radar signal received within a 
particular range cell index n given below (Chen [2]). 

+ (6.4) 

k=\ V C j 

where Mn is the total number of scattering points that share the same range, the 
target aspect angle 6(t) is expressed similarly to (2.3), and R(t) is given by (5.8) 
as shown below. 

6{t) = 6Q+Q}t + ^yt^ +.... (6.5) 


R{t) = R,+V,- " ^ \P{t\R^)dt' 

t 


( 6 . 6 ) 


where Ro is the range at the start of data collection (integration period) and Vt is 
the initial radial velocity. From Chapter II, the cross-range cell migration is 
maximal when the initial angle is zero. Setting 0() = 0 and assuming that there is 
no angular acceleration yields 6(t) = cot. Substitute (6.6) into (6.4) and, for zero 
initial velocity (as in our demonstrated case), we get 


;(0L =£exp< 




k=\ 


Ro- 




\P{t\R^)dt ' -I- xcos cot - yy. sin cot 


(6.7) 


By expanding sinter = (yr-(yy /6-h...and co^cot = \-co^t^ 12+... in (6.7) and 
ignoring the higher order terms, we obtain 


k=\ y c 

M 


A„ sin c 

(R, + x)-yyCDt- \P(t\RJdt' 


M. 


G t 


( 6 . 8 ) 


k=\ 
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It is easier to quantify the effect on cross-range cell migration if we 
consider the case where there is only one scattering point. The Doppler shift in 
the range cell with index n then can be obtained by differentiating the range 
related phase terms in (6.8) to give 


where 


at 


dt c 


A„sm6„ r 

(Ro + x)-y^a)t —-1 P(t)dt ' 


M. 


G t 


2/. 


A„ sin 6„ 

Mr- 


(6.9) 



for t<t^ 

for <t <t^+t^ (6.10) 


for t,+t^<t 


and po, td, ta and a depend on Rc, the distance between the blast point and the 
center of mass of the aircraft. Note that the second term Ap sin Ob P(t, Rc) / Mg in 
equation (6.9) represents the impulse response of the aircraft, and since the blast 
duration is much shorter than the integration period, this term is independent of 
T. Using Equation (2.40) with Arc = c / 2fcCoT, we can estimate the amount of 
cross-range cell shift to be 


AM = abs 


= abs 


■)-{U0)-UT)} 

1 [ 2f^ A^sinO, ^^ 

1 c M^ 



A„ sin 6„ 
- -HA 

coMq 


( 6 . 11 ) 
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where I/A can be obtained using (3.8) for known po, td and a. From the simulation 
parameters given in the program in Appendix A — po = 556274.3N/m\ td = 
O.OMlsec, a = 2.53, we obtain I/A = I974.7N-sec/rn. 

Also, in the first demonstration case, we have Arc = 1.23m, Ap = 64m\ Ob = 
-90, CO = 3.77deg/sec = 0.0658rad/sec. Mg = 22000kg. Using Equation (6.11), we get 
AM = 72 cells. This is more than the total number of cells in the cross-range. 
Hence, the scattering points will crossover into the subsequent time-history 
series of the image and ‘appear’ as if they belong the scattering points in that 
time-history. In this case, AM - M = 72 - 64 = 8 cells, which is about the same 
number of cells that was observed in the first demonstrated case. From Equation 
(6.11), the extent of the cross-range cell migration seems to depend on the 
aircraft impulse response to the blast wave. This has been verified also for other 
blast orientation and lower blast magnitude. 
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VII. CONCLUSION AND RECOMMENDATIONS 


A. CONCLUSION 

The work presented in this thesis is a preliminary demonstration of the 
possibility of using ISAR to conduct battle damage assessment (BDA) (based on 
computer simulation). The significant advantage of ISAR — as opposed to other 
types of sensors — is its ability to produce high resolution target images at 
greater surveillance ranges than optical systems. With high quality images, the 
positive evaluation of the extent of the target damage can be quickly made and 
immediate follow-on action can be taken if required. 

However, most ISAR images are generally static in nature and require 
relatively long processing periods to construct. The quality of the image is often 
tied to the processing time and the prior conditions associated with the 
processing. The foremost of these conditions is the assumption that the received 
radar signal is reflected from a target with only constant angular velocity 
components (an assumption that is used to map the Doppler variation of the 
respective scattering points caused by the angular rotation to the cross-range 
position of the scattering points). The presence of the radial velocities and 
deviation of the angular rotation rate will lead to image distortion. Unfortunately, 
these variations cannot be ignored during BDA. 

In this study, our main interest is on aircraft. When an aircraft is hit by 
explosive munitions two distinct behaviors can be observed. First, the blast wave 
of the explosion will interact with the surface body of the aircraft and cause the 
aircraft to have translational motion, rotational motion and vibration responses. 
The characteristics of these responses depend on the magnitude of the blast 
loading, the orientation of the aircraft surface to the blast wave and the direction 
of these responses with respect to the ISAR observing the target. Second, most 
explosive munitions against airborne targets carry pre-fabricated fragments and 
shrapnel aimed to penetrate the thin layer of the target to cause a variety of 
damages: from loss of flight control, tearing of structures to fuel tanks and 
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munitions detonation. The nature of this latter effect is more difficult to account 
for in simulation. 

1. Approach to Model the Aircraft Response When Hit 

Our approach was to examine the response of the aircraft in a 
deterministic manner by looking at blast loading effect. Even so, there was a 
need to narrow down the scope of the study to within the time frame of this thesis 
work and our main focus was solely on the translational and rotational motions of 
the aircraft. This is a reasonable approach if we consider the aircraft as a rigid 
body in air. These motions are likely to be the dominant responses apart for the 
latter effect caused by the fragment penetrations. Besides, they are directly 
detectable by the ISAR system. 

Our simulation begins with the formation of a blast wave in air and 
requires its characteristic parameters to be understood. They are discussed in 
the initial sections of Chapter III. It was shown that the blast wave can be 
represented by an impulse pressure force as a function of time by Equation (3.1) 
that depends on four key parameters — its peak overpressure, blast duration, 
arrival time and its waveform parameter. These four parameters are, in turn, 
dependent on the distance between the center of blast and the point of 
interaction between the pressure force and the surface. The parameters are 
usually described using empirical equations. 

The overall translational and rotational responses of the aircraft from the 
blast loading can be derived using the equation of motions (Newton’s 2"^ law). 
However, the non-symmetrical geometry of the aircraft surface and the unequal 
distribution of the blast wave (as a result of its spherical divergence) increase the 
complexity of the problem and it has been shown in Chapter III that it is difficult to 
yield an analytical closed-form solution. However, the results of Equations (3.21) 
and (3.22) can be obtained through the use of extensive numerical simulation. 
The methodological approach for the simulation was also briefly outlined in the 
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chapter. Due to the time constraint of this thesis work, we were (regrettably) not 
able to realize it. 

Instead, we simplified the study by modeling the translational response of 
the aircraft alone. This restriction can be accomplished with reasonable 
assumptions made concerning the nature of the blast wave and the surface area. 
We have assumed a planar blast wave with equal pressure distribution across 
the aircraft body. Also, the surface area of the aircraft during the evaluation 
period does is assumed to not change with time. These assumptions help to 
decouple the time and geometrical surface dependent pressure distribution of 
(3.21) into two separable variables shown in Chapter IV and in (5.1) — both 
independent of each other. The translational motion can be obtained in terms of 
radial acceleration, velocity and its resultant range variation with respect to time, 
by solving (5.1) analytically. However, it is not possible to similarly obtain the 
angular velocity variation since the resultant changes in angular momentum 
depend on the differences in pressure distribution across the surface. 

2. Approach to Overcome the Image Distortion 

The blast loading on the aircraft induces variations in both the radial and 
angular velocities. The velocity variations are encoded as changes in the 
Doppler frequencies of each returned scattering signal. This information is 
captured by the ISAR and use for radar imaging. However, the time resolution 
(which is equivalent to its integration time) of the conventional Fourier method of 
image processing is extremely poor and it is difficult to extract the dynamic 
frequency changes during and after the blast event occurs. As a result, the 
image becomes distorted due to the inability of the Fourier processing to account 
for the velocity variations. 

On the other hand, the Time-frequency Transform (TFT) method of signal 
analysis offers an added dimension allowing the signal to be represented as its 
frequency changes with time. We have shown through the use of a simple TFT 
method of analysis (which is the Short-time Frequency Transform or STFT) that 
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the dynamic behavior of the aircraft as it is being hit can be revealed in the ISAR 
images. The TFT images provide more information about the event than with 
FFT processing alone. The use of TFT for ISAR imaging is also beneficial to 
similar future work because of its ability to extract time-varying frequencies since 
many of these efforts will, in some way or another, involve examining Doppler 
variations of measured radar signals. 


B. RECOMMENDATIONS FOR FURTHER WORK 

Although, the use of ISAR systems with TFT processing has been 
demonstrated in this thesis for BDA, the simulation model used for proof-of- 
concept gives only an initial indication of this possibility. Further work is highly 
recommended to refine the present study in the three major areas discussed 
below. 


1. Complete Model of the Blast Loading on the Aircraft 

In this simulation, only the simplest case — the translational motion of the 
aircraft — is accounted for. Even then, there were certain assumptions made to 
reduce the magnitude of the problem so that some rough estimates could be 
obtained (for example, the possible angular variations were not considered). To 
enable a complete description of the aircraft response to the blast, there is a 
need to develop a full model that accounts for all the parameters given in 
Equations (3.21) and (3.22). It is not an impossible task, but requires careful 
consideration of the geometry of the problem between the aircraft (its surface 
distribution and orientation), the location of the blast point and its pressure 
distribution in time and space, and the position of the observing radar. The 
problem can also be extended from the present two-dimensional analysis to a 
real world three-dimensional case (but also with an increasing level of 
complexity). As discussed in Chapter III, this blast effect model can be developed 
as a separate model and its results can be simply expressed as a collection of 
time varying range and angular information over the entire integration period 

(following the methodological approach proposed). 
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2. Study of Target Break Up Process 

A further extension of the blast model is to include the target break up 
process due to the penetration of the fragments into the aircraft. However, the 
break up process is generally treated statistically and depends on the details of 
the fragment penetration. As a crude start, we could divide the aircraft into 
different sections and each section could be driven by different translational and 
rotational responses induced by the blast effect due to the difference in pressure 
wave distribution in each section. In this case, a series of range and angular 
variations for different sections can serve as input to the ISAR imaging model. 
The ISAR imaging model will then introduce the appropriate phase variations 
according to the range and angular inputs to respective scatterers categorized by 
the different sections. 

3. Improving the Time and Frequency Resoiution of TFT 

While the linear Short-time Frequency Transform (STFT) is the simplest to 
implement for our demonstration, we have seen that the time and frequency 
resolutions still lack the sensitivity required to provide faithful tracking of the 
aircraft dynamic behavior due to blast and simultaneously provide high quality 
radar imagery. Depending on the choice of the sliding window width and shape, 
the time-resolution can be adjusted in such a way that the dynamic motion of the 
target within the integration period can be extracted and tracked. But a major 
drawback of this approach is that, by gaining time resolution, the frequency 
resolution will necessarily suffer — and the quality of the image depends on the 
frequency resolution. This limitation is fundamental and is imposed by the 
uncertainty principle of time and frequency in any Fourier process. For a linear 
TFT, the best resolution is achievable using Gaussian window with the time- 
bandwidth product of y 2 . 

From a system point of view, one can further improve both resolutions by 
having a greater number of range and cross-range cells to represent the image, 

although other radar considerations (such as the carrier frequency, pulse interval 

91 



and frequency step size with respect to the size of the target and its rotation rate) 
will have to be factored in as they are discussed in Chapter II. 

The alternate approach is to work on enhancing TFT image processing. 
Most efforts involving time-frequency analysis utilize the more complex variations 
of STFT to overcome the time-frequency resolution limit imposed by the STFT. 
Examples of some of these TFT schemes include the Continuous Wavelet 
Transform (CWT), the bilinear Wigner-Ville Distribution (WVD), the Cohen’s class 
series and adaptive time-frequency representations. Many of these TFT 
approaches have reported better time-bandwidth products and, of these, the 
WVD seems to give best time-frequency resolution (but at the expense of 
significant cross interference terms). 
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APPENDIX A: MATLAB CODES 


% TFT ISAR IMAGE SIMULATOR 

% This program is developed to demonstrate the possibility of reconstructing 
% ISAR image using Time-Frequency Transform Method. It can be used for 
% viewing the dynamic behaviour of the target due to blast effect. 

% It contains 4 parts: 

% a. Part I - Define the radar parameters and Target Model. 

% b. Part 2 - Generate the blast effect on traslational motion. 

% c. Part 3 - Creates the returned stepped frequency signal. 

% (This is modified from the original code contained in Jae S. S. et 

% al "Range-Doppler Radar Imaging and Motion Compensation", Norwood, MA, 

% Artech House, 2001). 

% b. Part 4 - Construct the ISAR Image Using 2D EFT and TFT. 

% Assumptions Made: Range aligned. Motion compensated. Zero receiving system delay 
% and constant amplitude scatterers (unity). 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% Clear Memory 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clear all; fig_con = 0; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% PART 1: Define the radar parameters and Target Model. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% Define radar parameters 


M = 64; 

N = 64; 
step_f = 2.2e6; 
fc = 3e9; 

Ro = 22567.2; 
Vro = 0; 
Acc_ro = 0; 
Teta_o = 0; 
wo = 3.77; 
alo = -0.03; 
c = 3e8; 

T2 = 150e-6; 
T1 = l/step_f; 


% Total number of bursts 
% Total number of pulses 
% Frequency step (Hz) 

% Starting frequency (Hz) 

% Initial range (meters) 

% Initial translational velocity (m/ses) 

% Initial translational aceleration (m/sec-sq) 
% Initial angle of rotation (deg) 

% Initial angular velocity (deg/sec) 

% Initial angular acceleration (deg/sec-sq) 

% Speed of EM wave (m/sec) 

% Pulse repetition interval (sec) 

% Pulse length (sec) 


wo = wo*pi/180; alo = alo*pi/180; Teta_o = Teta_o*pi/180; 

T1 = l/step_f; 

% Define two dimensional density reflectivity function of Aircraft 
% Rows are x and y coordinates, 

Tgt_Model = zeros(N,M); 

Tgt_Model(29,14)=1;, Tgt_Model(31,14)=1; 

Tgt_Model(27,16)=!;, Tgt_Model(3 3,16)= 1; 

Tgt_Model(27,19)=!;, Tgt_Model(3 3,19)= 1; 

Tgt_Model(27,23)=l;, Tgt_Model(33,23)=l; 

Tgt_Model(25,25)=l;, Tgt_Model(35,25)=l; 
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Tgt_Model(23,27)=l 
Tgt_Model(21,29)=l 
Tgt_Model(19,31)=l 
Tgt_Model(17,33)=l 
Tgt_Model(15,35)=l 
Tgt_Model(13,37)=l 
Tgt_Model(15,38)=l 
Tgt_Model(17,39)=l 
Tgt_Model(19,38)=l 
Tgt_Model(21,37)=l 
Tgt_Model(23,36)=l 
Tgt_Model(25,35)=l 
Tgt_Model(27,34)=l 
Tgt_Model(27,37)=l 
Tgt_Model(27,40)=l 
Tgt_Model(28,43)=l 
Tgt_Model(26,45)=l 
Tgt_Model(24,47)=l 
Tgt_Model(28,47)=l 
Tgt_Model(30,42)=l 


, Tgt_Model(37,27)=l 
, Tgt_Model(39,29)=l 
, Tgt_Model(41,31)=l 
, Tgt_Model(43,33)=l 
, Tgt_Model(45,35)=l 
, Tgt_Model(47,37)=l 
, Tgt_Model(45,38)=l 
, Tgt_Model(43,39)=l 
, Tgt_Model(41,38)=l 
, Tgt_Model(39,37)=l 
, Tgt_Model(37,36)=l 
, Tgt_Model(35,35)=l 
, Tgt_Model(33,34)=l 
, Tgt_Model(33,37)=l 
, Tgt_Model(33,40)=l 
, Tgt_Model(32,43)=l 
, Tgt_Model(34,45)=l 
, Tgt_Model(36,47)=l 
, Tgt_Model(32,47)=l 


fig_con = figure(fig_con+l); 
imagesc(-Tgt_Model), colormap(gray) 
ylabel('Range Bin'), xlabel('Cross-Range Bin') 

titie('Density Reflectivity Distribution of the Target Scattering Points') 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% Part 2: Generate the blast effect on traslational motion. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


% Aircraft Area Profile - Front and Side 
Areafront = 7.6416; % Front Area in m-sq 
Areaside = 64; % Side Area in m-sq 
Mass_Aircraft = 22000; % Mass of Aircraft in kg 

% Blast Parameters - Air Explosive Profile. 

% Based on Medium Range SAM - Explosive Weight = 35kg, Miss Distance = 5m, 

% Target Height = 3000m, TNT Equivalent =1.5 (HE) 
fract = 0.25; 

Ar_Delay = fract*N*M*T2; % Arbitary Delay in fraction of integration time 
% Obtain from Excel Table - "Blast Effect Study - Input Parameters.xls" 

Po = 556274.3; % Peak Overpressure in Pa (N/m-sq) 

ta = 0.015593 h- Ar_Delay; %arrival time in sec 
td = 0.014116; %blast duration in sec 
alpha = 2.53; %waveform parameter 

% Geometry of Blast - Equation (5.10), (5.7), (5.8) & (5.9) 

TetaB_deg = -90; % direction of blast (deg) CW ref to neg. crossrange axis 

TetaB = TetaB_deg / 180 * pi; 

AreaP = abs(Areafront * cos(TetaB H- Teta_o)) H- abs(Areaside * sin(TetaB H- Teta_o)); % Projected Area 
Eyo = Po * AreaP * sin(TetaB); 

Ayo = - Eyo / Mass_Aircraft; 

% Sampling times Si 
% Sik vector of 64x64 elements 
i= 1;M*N; 
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Sik = (Tl/2) + 2*Ro/c + ((i-l)*T2); % Using Equation (2.9) 


% Range Variation ('!' due to blast effect, '0' linear expression) 

With_Blast = 1; % Blast Switch 

if With_Blast % Using Equation (5.7), (5.8) and (5.9) 

Rik = zeros(l,N*M); 
for vect_index = 1 ;M*N 
if (Sik(vect_index) <= ta) 

Rik(vect_index) = Ro + Vro*Sik(vect_index); 
elseif (Sik(vect_index) > ta) && (Sik(vect_index) <= (ta + td)) 

Rik(vect_index) = Ro + Vro*Sik(vect_index) + Ayo*((td/alpha)^2 * (2/alpha-l)*(l-exp( 
-alpha*(Sik(vect_index)-ta)/td)) + td/alpha''2 * (Sik(vect_index)-ta)*( 
alpha-1- exp(-alpha*(Sik(vect_index)-ta)/ td))); 

else 

Rik(vect_index) = Ro H- Vro*Sik(vect_index) + Ayo*((td/alpha - td/alpha^2 * (1 - exp( 
-alpha)))*(Sik(vect_index) - (ta + td)) + td^2/alpha^3 * (2 - alpha) * ( 

1 - exp(-alpha)) + td''2/alpha^2 * (alpha - 1 - exp(-alpha))); 

end 

end 

else 

Rik = Ro H- Vro*Sik + 0.5*Acc_ro*Sik.''2; % Equation (2.2) 

end 

Rik = reshape(Rik,N,M); 

% Angular positions - Equation (2.3) 

Tetik = Teta_o + wo*Sik + 0.5*alo*Sik.^2; 

Tetik = reshape(Tetik,N,M); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% Part 3: Creates the returned stepped frequency signal. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Erequency steps - Equation (2.8) 
i= 1:M; 

fil = fc H- (i-l)*step_f; fi = []; 
for i = 1:N fi = [fi fil]; end 
clear fi 1; 

fis = reshape(fi,N,M); 

% Create the target image by summing all the returns together. 

% Step 1: Eind the number of point scatterers on the image 
numpo = length(find(Tgt_Model == 1)); 

[xx yy] = find(Tgt_Model == 1); 
pos = 1; 

target = zeros(M,N); 

% Step 2; Perform the summation - Equation (2.10) 
for i = 1: numpo 

target = target + exp(-i*4*pi/c*fis.*Rik).*exp(-i*2*pi*( (xx(pos)-N/2)*2*(fis/c).*cos(Tetik) 

- (yy(pos)-N/2)*2*(fis/c).*sin(Tetik))); 
pos = pos H-1; 
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% PART 4: Construct the ISAR Image Using 2D PET and TET. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

% TET Cube Initialisation 

stft_cube = zeros([N M N]); % Instant Prequency-Time History-Range Cell 
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% Obtain the Synthetic Range Profile of the target - Equation (2.11) 
th_target = fft(target); % fft on column of matrix 

% Construct the Conventional FFT Image - Equation (2.12) 

fig_con = figure(fig_con+l); 

imagesc(fftshift(abs(fft(th_target.').'))); 

ylabel('Range Bin') 

xlabel('Cross-Range Bin') 

colormap(gray); 

% Perform TFT on the doppler history of each Range cell index using STFT 
% Ref: Chapter IV Section D 
% Window function 

h = window(45,'gauss',0.005); % For gauss - need to adjust width according to sigma / 0.005 
for vect_index = 1 :M 

stft_cube(:,:,vect_index) = tfrstft(th_target(vect_index,:).',l:M,M,h).'; 
end 

% Generate and Plot the Reconstructed ISAR 
for vect_indexl = 1:4 

fig_con = figure(fig_con+l); 
for vect_index =1:16 

fii = fftshift(squeeze(stft_cube(vect_index+(vect_indexl-l)*16,:,:)).'); 
subplot(4,4,vect_index), imagesc(abs(fii)), colormap(gray); 

THmarker = vect_index+(vect_indexl-l)*16; 

title(['TH',num2str(THmarker),' (t = ', num2str(THmarker*T2*N),'sec)']); 
axis off; 
end 
end 
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